OPTIMAL  PERIODIC  INSPECTION  OF  A  STOCHASTICALLY 

DEGRADING  SYSTEM 

THESIS 

Timothy  B.  Booher,  Captain,  USAF 
AFIT  /  GOR/ENS/ 06-04 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Patterson  Air  Force  Base,  Ohio 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the 
official  policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense  or 
the  United  States  Government. 


AFIT  /  GOR/ENS  /  06-04 


OPTIMAL  PERIODIC  INSPECTION  OF  A  STOCHASTICALLY 

DEGRADING  SYSTEM 


THESIS 


Presented  to  the  Faculty 
Department  of  Operational  Sciences 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
in  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science  in  Operations  Research 


Timothy  B.  Booher,  S.B. 
Captain,  USAF 


March  2006 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


AFIT/GOR/ENS/06-04 


OPTIMAL  PERIODIC  INSPECTION  OF  A  STOCHASTICALLY 

DEGRADING  SYSTEM 

Timothy  B.  Booher,  S.B. 

Captain,  USAF 


Approved: 


i  6i 

Date 

L  m.  °L 


Lt  Col  Mark  Abramson,  Ph.D. 
Committee  Member 


Date 


AFIT  /  GOR/ENS  /  06-04 


Abstract 

This  thesis  develops  and  analyzes  a  procedure  to  determine  the  optimal  inspect¬ 
ion  interval  that  maximizes  the  limiting  average  availability  of  a  stochastically  de¬ 
grading  component  operating  in  a  randomly  evolving  environment.  The  component 
is  inspected  periodically,  and  if  the  total  observed  cumulative  degradation  exceeds 
a  fixed  threshold  value,  the  component  is  instantly  replaced  with  a  new,  statisti¬ 
cally  identical  component.  Degradation  is  due  to  a  combination  of  continuous  wear 
caused  by  the  component’s  random  operating  environment,  as  well  as  damage  due 
to  randomly  occurring  shocks  of  random  magnitude.  In  order  to  compute  an  opti¬ 
mal  inspection  interval  and  corresponding  limiting  average  availability,  a  nonlinear 
program  is  formulated  and  solved  using  a  direct  search  algorithm  in  conjunction 
with  numerical  Laplace  transform  inversion.  Techniques  are  developed  to  signifi¬ 
cantly  decrease  the  time  required  to  compute  the  approximate  optimal  solutions. 
The  mathematical  programming  formulation  and  solution  techniques  are  illustrated 
through  a  series  of  increasingly  complex  example  problems. 
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OPTIMAL  PERIODIC  INSPECTION  OF  A  STOCHASTICALLY 

DEGRADING  SYSTEM 

1.  Introduction 

All  systems  experience  gradual  deterioration  over  time  until  they  must  ulti¬ 
mately  be  repaired  or  replaced.  In  order  to  maintain  these  systems,  engineers  and 
analysts  must  decide  when  and  how  to  inspect  and  maintain  them.  The  study  of 
optimal  maintenance  planning  provides  a  quantitative  basis  for  maintenance  deci¬ 
sions  and  is  of  interest  to  industrial,  governmental,  and  military  organizations.  In 
most  typical  scenarios,  optimal  maintenance  is  concerned  with  maintaining  a  system 
in  a  manner  that  maximizes  some  measure  or  benefit,  or  minimizes  the  long-run 
average  cost.  Models  for  optimal  maintenance  planning  often  employ  probability 
and  stochastic  process  theory  in  an  attempt  to  model  realistic  complexities  which 
are  inherent  in  many  components,  sub-systems,  and  systems. 

Until  recently,  many  optimal  maintenance  models  have  lacked  sufficient  detail 
to  account  for  the  complex  interactions  of  degradation  mechanisms  that  determine 
component  lifetime  distributions.  Consequently,  many  of  the  potential  gains  of  op¬ 
timal  maintenance  models  remain  untapped.  While  the  models  considered  in  this 
thesis  are  general,  and  therefore  apply  to  a  wide  array  of  scenarios,  the  application 
area  of  focus  for  this  research  is  the  United  States  Air  Force’s  maintenance  policies. 

The  United  States  Air  Force  maintains  equipment  ranging  from  multi-billion 
dollar  weapons  platforms,  such  as  the  B-2  stealth  bomber,  to  the  most  mundane  tools 
and  facilities.  Optimizing  the  inspection  intervals  for  Air  Force  equipment  is  a  fertile 
area  with  the  potential  to  reduce  costs  dramatically  while  simultaneously  increasing 
operational  availability.  Optimal  maintenance  theory  is  especially  needed  in  the 
current  environment  where  senior  Air  Force  leadership  has  challenged  all  Air  Force 
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maintenance  processes  to  increase  availability  by  an  average  of  20%  while  achieving 
a  simultaneous  10%  decrease  in  cost  by  the  year  2011.  This  research  provides  a 
previously  non-existent  capability  to  link  availability  and  cost  together  in  a  single 
model  producing  an  approximate  optimum  inter-inspection  interval  that  maximizes 
the  limiting  average  availability  while  satisfying  a  prespecified  budget. 

Current  U.S.  Air  Force  policies  mandate  the  implementation  of  reliability- 
centered  maintenance ,  which  specifies  that  all  Air  Force  organizations  maintaining 
equipment  must  have  a  plan  to  consider  the  reliability  measures  of  the  components 
for  which  they  are  responsible  to  maintain.  Most  often,  reliability  engineers  in  these 
organizations  use  statistical  data,  often  collected  from  an  analogous  system,  to  fit 
a  Weibull  distribution  to  estimate  the  necessary  reliability  measures  which  are  then 
used  to  determine  inter-inspection  intervals.  The  cost  of  this  maintenance  policy 
is  determined  at  some  later  time,  and  an  iterative  refinement  process  ensues  before 
converging  on  a  compromise  between  cost  and  availability.  Moreover,  the  current 
methods  used  to  determine  if  existing  inspection-interval  lengths  can  be  extended  are 
even  more  simplistic.  An  informal  survey  found  that  the  most  common  method  used 
by  Air  Force  depot  engineers  to  extend  inspection  intervals  is  to  set  aside  several 
systems  for  inspection  at  the  increased  interval  of  interest.  When  inspected  after 
this  prolonged  period  of  time,  the  systems  were  closely  inspected  to  see  a  posteriori 
if  any  damage  was  outside  acceptable  limits. 

For  currently  fielded  systems,  system  managers  acknowledge  that  90%  of  the 
reliability  characteristics  of  a  system  are  set  by  the  time  10%  of  the  program  dollars 
are  expended.  Therefore,  it  is  critical  that  analytical  tools  are  available  as  early 
as  possible  in  the  acquisition  process.  While  the  assumptions  behind  any  model 
(mathematical  or  otherwise)  cannot  be  completely  validated,  an  educated  initial 
value  for  an  appropriate  inter-inspection  interval  that  maximizes  availability  would 
be  a  useful  starting  point  for  engineers.  An  application  area  of  particular  interest 
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to  the  Air  Force  is  to  use  this  starting  point  when  determining  inspection  intervals 
during  the  logistics  supportability  analysis  required  before  a  system  can  be  fielded. 

An  accurate  understanding  of  a  component’s  reliability  has  special  meaning 
in  a  military  context.  Military  applications  require  equipment  to  operate  in  the 
harshest  conditions  and  must  have  reliability  measures  sufficient  to  both  protect  lives 
and  accomplish  the  required  mission.  In  order  to  illustrate  and  further  motivate 
the  relevance  of  this  thesis,  an  application  concerning  the  maintenance  of  Radar 
Absorbing  Material  (RAM)  is  next  presented. 

Low  observable  (LO)  technologies  make  weapons  systems  difficult  to  detect, 
track,  and  engage.  These  weapons  systems,  termed  stealth  assets,  are  vital  to  the 
future  of  the  Air  Force  and  have  been  combat-tested  in  several  operations  includ¬ 
ing  Operation  Desert  Storm,  Operation  Allied  Force  and  Operation  Iraqi  Freedom. 
Stealth  aircraft  are  particularly  valuable  against  high-value  targets,  which  are  often 
either  out  of  range  of  conventional  aircraft  or  too  heavily  defended  for  conventional 
aircraft  to  strike.  This  critical  role  in  combat  operations  will  increase  in  impor¬ 
tance  in  the  coming  years,  as  Russian-built  surface-to-air  missile  systems,  such  as 
the  SA-10,  SA-12,  SA-15,  and  SA-20  appear  on  the  open  market.  These  integrated 
air  defense  systems  are  highly  mobile,  networked,  and  possess  an  entire  suite  of  anti¬ 
aircraft  engagement  mechanisms.  Because  of  these  capabilities,  surface-to-air  missile 
strikes,  and  counter-ground  jamming,  combined  with  conventional  strike  aircraft, 
will  not  be  adequate  to  gain  and  maintain  air  supremacy.  Clearly,  low  observable 
technology  will  play  a  vital  role  in  both  low  and  high-intensity  operations. 

While  LO  technology  mitigates  detection  from  many  sensors,  such  as  heat 
seekers  (infrared),  sound  detectors,  and  even  the  human  eye,  the  ability  to  reduce  the 
radar  cross  section  (RCS)  is  a  key  component  of  stealth.  There  are  two  approaches 
to  reduce  the  passive  radar  cross  section:  shaping  to  minimize  backscatter,  and  using 
RAM  coating  for  energy  absorption  and  cancelation.  While  both  these  mechanisms 
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have  proven  effective  in  reducing  RCS,  the  maintenance  requirement  associated  with 
RAM  is  often  an  area  of  concern. 

Not  only  are  these  materials  extremely  sensitive  to  the  ambient  environment, 
they  must  maintain  precise  tolerances  for  an  entire  array  of  electrical  properties  in 
order  for  the  entire  system  to  effectively  reduce  a  system’s  RCS.  Therefore,  the  effect 
of  environmental  exposure  has  a  direct  impact  on  aircraft  RCS  and  consequently  air¬ 
craft  survivability.  This  concern  persists  even  in  the  most  benign  environments,  as 
optimal  flying  conditions  have  been  shown  to  induce  noticeable  wear  on  RAM  [137]. 
Sunny  weather,  rain,  and  sandy  environments  all  cause  RAM  to  degrade  linearly  at 
different  rates.  The  B-2,  for  example,  must  be  stored  in  a  climate-controlled  hangar 
to  mitigate  RAM  near-field  reflectivity  loss.  Near-field  reflectivity  is  an  important 
measure  of  RAM  performance  and  can  be  measured  at  a  specific  location.  Measur¬ 
ing  near-field  reflectivity  per  periodic  inspections  by  maintenance  personnel  is  more 
cost  effective  than  measuring  the  overall  RCS  degradation,  which  requires  an  entire 
aircraft  to  be  tested  in  a  specialized  facility. 

Environment-dependent  wear  alone  cannot  adequately  characterize  the  RAM 
deterioration  process.  Shocks  also  provide  an  extremely  important  degradation 
mechanism.  For  example,  in-air  refueling  stops  are  short  periods  of  time  that  can 
cause  rapid  damage  when  the  refueling  aircraft’s  boom  scratches  the  sensitive  LO 
surfaces.  Other  shocks  include  bird  strikes  and  maintenance  personnel  touching  the 
RAM  surface.  The  linear  wear  rates  for  near-field  reflectivity  loss  are  easily  mea¬ 
sured  experimentally.  Shock  damage  is  also  measurable  from  individual  maintenance 
records,  and  shock  arrival  rates  are  readily  obtained  from  historical  data  in  current 
maintenance  databases. 

Due  to  the  deterioration  process  described  above,  RAM  must  be  routinely 
inspected  and  replaced.  Since  RAM  deterioration  can  only  be  measured  when  in¬ 
spected,  failures  are  said  to  be  hidden  or  non  self-announcing.  The  inspection  inter¬ 
val  is  very  sensitive  with  regard  to  cost  and  operational  availability  since  in  order  to 
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inspect  RAM  deterioration,  an  aircraft  must  be  removed  from  operational  use.  An 
overly  small  inspection  interval  could  be  operationally  and  financially  expensive,  but 
too  lax  an  inspection  interval  could  potentially  result  in  a  compromise  of  aircraft 
survivability  and  risk  of  loss  of  an  extremely  costly  asset. 

A  stochastic  model  developed  by  Kharoufeh  et  al.  [67]  is  well-suited  for  the 
RAM  degradation  process  and  presents  reliability  and  availability  measures  for  a 
component  subject  to  degradation  from  a  simultaneous  combination  of  linear  wear 
from  a  random  environment  and  from  random  shocks.  The  combination  of  these 
two  degradation  mechanisms  gives  this  model  considerable  flexibility  in  addressing  a 
wide-range  of  applications  as  many  systems  are  maintained  by  inspections  conducted 
according  to  a  fixed,  deterministic  interval.  Concerning  the  previous  example,  with 
correctly  identified  parameters  for  the  shock  and  wear,  and  an  appropriate  opti¬ 
mization  model,  it  is  possible  to  compute  an  optimal  inter-inspection  interval  that 
maximizes  the  limiting  average  availability  of  a  stealth  weapons  system. 

Availability  is,  in  the  most  general  sense,  the  proportion  of  time  that  a  compo¬ 
nent  or  system  is  available  to  be  used  for  its  intended  purpose.  Naturally,  availability 
is  an  important  measure  of  concern  in  many  applications  and  is  generally  desired  to 
be  as  high  as  possible.  While  the  model  of  Kharoufeh  et  al.  [67]  has  wide  applicability 
and  provides  a  mechanism  to  understand  the  impact  of  a  particular  selection  of  the 
inter-inspection  interval  on  reliability  and  availability  measures,  no  previous  research 
has  developed  a  methodology  to  maximize  the  limiting  average  availability  through 
selection  of  the  inter-inspection  interval  for  a  model  that  considers  degradation  due 
to  environment- induced  wear  and  random  shocks. 

1.1  Problem  Definition  and  Methodology 

This  thesis  will  develop  an  appropriate  cost  function  and  formal  nonlinear 
program  to  maximize  the  limiting  average  availability  of  a  component  subject  to 
linear  wear  and  random  shocks.  The  decision  variable  of  interest  is  the  fixed  length  of 
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the  inter-inspection  interval.  An  appropriate  optimization  technique  will  be  selected 
and  implemented  through  an  accessible  set  of  computer  codes.  Finally,  several  test 
cases  will  be  constructed  to  illustrate  the  optimization  technique  in  practice. 

However,  in  order  for  the  methodology  to  be  useful  in  an  industrial  or  mili¬ 
tary  setting,  it  must  be  able  to  obtain  an  answer  in  a  reasonable  amount  of  time. 
Therefore,  the  research  objectives  of  this  thesis  are:  (1)  to  construct  the  necessary 
optimization  methodology  to  maximize  availability,  and  (2)  to  produce  this  answer 
in  the  least  amount  of  time  possible. 

To  accomplish  these  two  objectives,  the  stochastic  degradation  model  devel¬ 
oped  by  Kharoufch  et  al.  [67]  will  be  reviewed.  Next,  a  cost  function  will  be  developed 
which  considers  the  long-run  cost  per  cycle  due  to  component  replacement,  inspec¬ 
tions,  and  downtime.  The  objective  function  of  [67]  and  the  associated  cost  function 
combine  to  form  a  nonlinear  program.  However,  because  the  representation  for  the 
limiting  average  availability  is  provided  only  in  the  form  of  a  Laplace  transform, 
real-domain  derivative  information  is  not  available,  and  a  solution  requires  the  use 
of  Laplace  transform  inversion  techniques.  For  this  reason,  many  standard  optimiza¬ 
tion  techniques  are  inadequate;  thus  derivative-free  algorithms  (particularly  pattern 
search)  are  considered. 

Since  the  objective  function  producing  the  limiting  average  availability  is  com¬ 
putationally  intensive,  a  study  was  conducted  to  reduce  the  run  times  of  its  most 
frequently  called  subfunctions:  the  matrix  exponential  and  the  Laplace  transform 
inversion  algorithm.  Four  methods  to  compute  the  matrix  exponential  were  studied 
and  a  Laplace  transform  inversion  algorithm  was  selected.  A  complexity  analysis 
was  conducted  on  the  direct  computation  of  the  limiting  average  availability  and 
cost  computations.  This  leads  to  the  development  of  a  method  that  dramatically 
reduces  run  times  as  the  number  of  environmental  states  and  system  maximum  life¬ 
time  increase.  All  computational  methods  and  the  overall  optimization  problem  are 
illustrated  through  five  test-cases. 
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1.2  Thesis  Outline 


The  next  chapter  includes  a  brief  survey  of  the  vast  literature  in  the  field  of  op¬ 
timal  maintenance.  Starting  with  the  first  transition  of  ideas  from  other  disciplines 
into  an  optimal  maintenance  context,  the  review  surveys  three  frequently  used  sto¬ 
chastic  deterioration  mechanisms  and  the  literature  on  optimal  replacement.  The 
literature  reviewed  will  describe  the  history  and  variations  of  currently  used  optimal 
maintenance  models,  with  an  emphasis  on  the  historical  thread  that  develops  into 
the  model  examined  in  this  thesis.  The  goal  of  this  survey  is  to  provide  the  reader 
the  necessary  background  to  understand  the  particular  contributions  of  the  thesis 
in  their  historical  context.  This  chapter  highlights  the  need  for  an  implementable 
optimization  methodology  to  determine  the  optimal  period  of  time  to  wait  between 
inspections  to  maximize  limiting  average  availability  for  a  given  budget  constraint 
and  associated  cost  function. 

The  formal  notation  and  mathematical  model  are  developed  in  chapter  3.  The 
first  section  describes  the  notation  and  formulation  of  the  stochastic  degradation 
model  developed  by  Kharoufeh  et  al.  [67],  which  presents  reliability  and  availability 
measures  for  a  system  that  degrades  according  to  a  simultaneous  combination  of 
linear  wear  and  random  shocks.  The  second  section  describes  the  development  of  a 
cost  function  that  considers  the  long-run  cost  per  cycle  of  each  inspection.  In  the 
third  section,  a  formal,  constrained  nonlinear  programming  formulation  is  presented. 

Chapter  4  discusses  the  solution  methodology  used  to  solve  the  nonlinear  pro¬ 
gram  presented  in  chapter  3.  This  chapter  starts  with  a  discussion  of  potential 
solution  methodologies  and  their  degree  of  applicability  to  the  nonlinear  program 
in  question.  In  particular,  non-derivative-based  search  procedures  are  considered 
with  pattern  search  chosen  as  the  best  overall  optimization  procedure.  In  order  to 
decrease  the  run  time  of  the  optimization  procedure,  algorithmic  refinements  are 
made  to  the  matrix  exponential  and  the  Laplace  transform  inversion  algorithm.  The 
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chapter  also  presents  the  development  of  a  numerical  method  to  dramatically  reduce 
run  times  for  the  chosen  methodology. 

Chapter  5  presents  the  results  of  several  numerical  experiments.  The  limiting 
average  availability  is  computed  using  pattern  search  methods  with  four  different 
methods  to  calculate  the  matrix  exponential.  The  result  from  this  experiment  is 
the  most  suitable  matrix  exponential  method,  which  is  then  used  in  a  second  ex¬ 
periment  to  isolate  the  gains  provided  by  a  new  method  developed  to  decrease  the 
computational  run  time.  The  thesis  is  concluded  in  chapter  6  by  summarizing  the 
main  results  and  discussing  the  contributions  of  this  research,  and  making  some  final 
recommendations  as  well  as  suggestions  for  future  research  directions. 
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2.  Review  of  the  Literature 

This  chapter  provides  an  overview  of  the  literature  related  to  the  area  of 
operations  research  termed  optimal  maintenance,  with  an  emphasis  on  optimal  re¬ 
placement  models  and  availability  optimization.  Optimal  maintenance  is  a  well- 
established  field  within  operations  research  concerned  with  maintaining  a  component 
in  a  manner  that  maximizes  profit  or  minimizes  cost. 

This  review  begins  with  a  discussion  of  motivations  for  implementing  optimal 
maintenance  and  its  general  methodology  followed  by  a  discussion  of  the  history  of 
optimal  maintenance  models.  Since  the  foundation  of  any  maintenance  model  relies 
on  the  underlying  deterioration  process  and  failure  behavior  of  the  component,  three 
common  degradation  models  are  reviewed:  shock  models,  wear  models  and  com¬ 
pound  models  that  consider  degradation  from  both  shock  and  wear.  After  covering 
degradation  models,  optimal  replacement  models  for  both  preventive  and  corrective 
maintenance  are  reviewed.  That  section  concludes  with  a  discussion  of  recent  models 
for  availability  analysis  of  complex  components. 

All  components  degrade  over  time  and  actions  must  be  taken  to  keep  equipment 
operating  for  its  intended  purpose.  These  actions  include  both  preventive  actions 
undertaken  at  regular  intervals  to  prevent  an  unacceptable  loss  of  performance  and 
corrective  actions  taken  to  restore  a  failed  component  to  an  operational  state.  In 
both  cases,  maintenance  actions  can  be  accomplished  at  fixed  intervals  or  based  on 
conditions  concerning  the  component  of  interest. 

In  this  review,  a  component  describes  any  system  or  subsystem  that  is  modeled 
and  evaluated  as  a  single  entity.  Thus,  a  system  that  consists  of  multiple  components, 
yet  works  towards  a  single  objective  can  be  considered  a  component  in  order  to 
reflect  the  aggregate  system-level  behavior.  Thus,  the  detailed  interactions  between 
the  various  subcomponents  can  be  neglected. 
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In  many  industries,  costs  associated  with  the  allocation  of  resources  to  mainte¬ 
nance  are  very  significant,  and  maintenance  personnel  comprise  a  significant  number 
of  the  total  workforce.  The  potential  impact  of  maintenance  at  the  level  of  opera¬ 
tions  and  logistics  is  considerable,  and  the  financial  implications  of  maintenance  can 
be  substantial.  Therefore,  models  are  required  to  gain  greater  insight  into  mainte¬ 
nance  processes  in  order  to  better  set  inspection  intervals  and  understand  process 
dynamics. 

Optimal  maintenance  models  have  diverse  application  areas,  including  crack 
growth  ([123]  and  [100]),  airframes  [112],  offshore  structures,  coastal  flood  barriers 
subject  to  erosion  [125],  and  many  other  industrial  and  military  settings.  While  ex¬ 
tensive  effort  has  been  devoted  to  developing  simulation  models,  as  surveyed  by  Cza- 
jkiewicz  [49],  simulation  modeling  does  not  give  insight  into  the  underlying  structure 
of  maintenance  problems  and  may  provide  inaccurate  answers,  while  simultaneously 
shielding  the  source  of  these  errors.  Further,  even  with  recent  advances  in  processor 
speed,  simulations  can  take  long  periods  of  time  to  compute  an  approximate  an¬ 
swer,  as  compared  to  the  seconds  required  to  compute  an  exact  answer  by  analytical 
methods. 

2.1  The  History  of  Optimal  Maintenance 

The  beginning  of  modern  optimal  maintenance  is  closely  correlated  to  the 
beginning  of  operations  research  in  general  which  was  developed  during  the  second 
World  War.  The  foundations  of  the  first  optimal  maintenance  models  came  from 
actuarial  research  performed  in  Switzerland  in  the  early  20th  century.  This  research 
focused  on  determining  the  number  of  annual  accessions  required  to  maintain  a  finite 
body  of  policyholders  [82],  Lotka  [82]  adapted  this  research  through  the  use  of  the 
theory  of  renewals  to  address  industrial  replacement  problems.  Lotka’s  research  was 
followed  by  more  applications  of  renewal  theory  and  actuarial  models  to  structural 
reliability  and  fatigue  failure  applications.  The  crossover  from  other  disciplines  into 
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optimal  maintenance  theory  continued,  as  theory  originally  developed  for  population 
analysis  and  problems  in  genetics  was  applied  to  a  number  of  industrial  problems. 

During  the  second  World  War,  the  fields  of  reliability  and  optimal  mainte¬ 
nance  were  studied  with  renewed  interest  as  both  the  industrial  economy  and  aca¬ 
demic  community  focused  on  assisting  the  war  effort.  During  this  time,  work  by 
Weibull  [142]  focused  on  approximating  probability  distributions  to  model  the  fail¬ 
ure  mechanics  of  materials  and  introduced  the  well-known  Weibull  distribution  for 
use  in  modeling  component  lifetimes.  Shortly  thereafter,  Davis  [50]  demonstrated 
that  the  exponential  distribution  worked  well  for  modeling  many  components.  Today, 
the  exponential  distribution  remains  the  distribution  of  choice  for  many  maintenance 
modeling  problems,  as  it  has  been  demonstrated  to  aptly  model  the  real-world.  It 
is  also  the  only  distribution  to  possess  the  memoryless  (or  Markov)  property,  which 
allows  for  easy  aggregation  of  failure  rates  of  subcomponents  to  determine  the  failure 
rate  of  the  overall  component. 

Maintenance  models  increased  dramatically  in  complexity  in  the  1960s  as  the 
study  of  preventive  maintenance  as  a  research  discipline  began  to  appear  in  the 
literature.  Researchers  such  as  Barlow  [28],  Proschan  [29],  Jorgenson  [65],  McCall 
[87],  Radner  [65],  and  Hunter  [28]  contributed  extensively  to  the  early  development 
of  maintenance  optimization  models.  Applications  in  this  era  focused  on  systems 
with  potentially  catastrophic  failures,  such  as  nuclear  power  plants  and  optimal 
maintenance  for  intercontinental  ballistic  missiles  [65].  The  increased  complexity  of 
many  of  these  models  required  the  relaxation  of  the  exponential  distribution,  and 
research  expanded  to  include  semi-markov  processes  [122], 

Current  research  uses  stochastic  processes  to  describe  the  failure-generating 
mechanisms  of  optimal  maintenance  models.  Renewal  theory  remains  the  prevalent 
method  used  to  model  stochastic  failure  processes.  Despite  the  broadening  scope 
of  modern  research  to  include  more  complex  components,  many  studies  have  shown 
that  renewal  theory  is  capable  of  accommodating  these  complexities.  The  following 
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sections  will  survey  the  current  literature  in  optimal  maintenance  starting  with  a 
review  of  the  underlying  deterioration  process  and  failure  behavior  of  components  in 
optimal  maintenance  models. 

2.2  Degradation  Mechanisms 

Stochastic  degradation  models  are  mathematical  models  which  attempt  to  de¬ 
scribe  a  component’s  deterioration  over  time.  There  are  two  primary  mechanisms 
presented  in  the  literature  for  modeling  a  component’s  aging  process  over  time: 
the  use  of  existing  lifetime  distributions  or  mathematical  modeling  of  the  physical 
dynamics  that  cause  the  failure  behavior.  Lai  and  Xie  [75]  provide  an  overview 
of  model  aging  and  dependence  characteristics  in  reliability  and  survival  analysis. 
Existing  distributions  fall  into  the  broad  categories  classified  by  how  their  failure 
rate  changes  over  time:  increasing,  decreasing,  or  constant.  Since  the  foundation  of 
any  maintenance  model  relies  on  the  underlying  degradation  process  of  the  system 
or  component  under  consideration,  the  following  section  will  review  three  primary 
degradation  models  actively  used  in  current  research:  shock  models,  wear  models 
and  compound  damage  models. 

Shock  models  are  effective  in  describing  a  component  whose  degradation  is  the 
result  of  a  collection  of  distinct  stresses,  termed  shocks,  applied  at  discrete  points  in 
time.  The  literature  provides  a  wide  variety  of  methods  to  describe  the  frequency  of 
occurrence,  and  magnitude  of  damage  caused  by  shocks.  Most  commonly,  however, 
probability  distributions  are  used  to  model  damage  magnitude,  and  shocks  arrive 
according  to  a  Poisson  process.  There  are  two  broad  areas  of  the  failure  mechanisms 
into  which  shock  models  can  be  classified:  cumulative  damage  shock  models  and 
maximum  shock  models.  Cumulative  damage  models  consider  a  failure  to  occur 
when  the  sum  of  the  effects  of  random  shocks  over  time  exceeds  a  particular  threshold 
value.  In  contrast,  components  described  by  maximum  shock  models  fail  when  the 
magnitude  of  a  single  shock  exceeds  a  particular  threshold  value. 
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Early  Poisson  shock  models  are  presented  by  Esary  et  al.  [54]  and  Ross  [110]. 
Shanthikumar  and  Sumita  [115]  present  a  more  complex  shock  model  in  which  shock 
magnitudes  are  correlated  with  the  length  of  the  interval  to  the  next  shock.  Two 
years  later,  Sumita  and  Shanthikumar  [129]  extended  their  previous  model  with 
the  development  of  a  cumulative  shock  model.  Rangan  and  Sarada  [108]  present  the 
earliest  paper  found  during  this  survey  that  uses  a  non-homogeneous  Poisson  process 
to  model  the  shock  arrival  process. 

Mallor  and  Santos  [84]  classify  general  shock  models  used  for  component  re¬ 
liability.  They  also  present  application  areas,  noting  that  extreme  and  cumulative 
shock  models  may  be  appropriate  descriptions  for  the  fracture  of  brittle  materials 
and  for  damage  due  to  earthquakes  or  volcanic  activity.  They  employ  Laplace  trans¬ 
forms  to  provide  the  distribution  function  of  the  component  failure  time  and  its 
mean  value.  Rade  [107]  presented  a  shock  model  with  a  finite  number  of  identical 
components,  each  receiving  shocks  arriving  occurring  to  a  Poisson  process.  He  used 
this  model  to  calculate  the  time-dependent  probability  to  failure  for  each  identical 
component.  Nakagawa  [96]  developed  a  replacement  policy  for  Rade’s  system  of 
identical  components,  which  exchanges  a  given  component  before  failure  if  the  total 
number  of  failed  components  is  more  than  a  fixed  number,  n,  and  which  replaces  the 
system  if  all  n  components  have  failed.  Nakagawa  determined  the  optimal  number, 
n,  that  minimizes  the  expected  cost.  Igaki  et  al.  [62]  consider  a  state-dependent 
shock  model  influenced  by  an  external  system,  in  which  both  the  interarrival  time 
and  the  magnitude  of  the  shock  are  determined  by  a  Markov  process.  Ebrahimi  [53] 
presents  a  model  subject  to  shocks  governed  by  a  Poisson  process,  where  damage 
accumulates  additively  and  the  component  fails  if  the  total  damage  exceeds  a  certain 
capacity  or  threshold.  His  results  allow  for  the  comparison  of  two  random  processes 
by  stochastic  ordering. 

While  shock  models  can  adequately  model  many  components,  they  require  the 
application  of  stresses  at  finite  intervals  of  time.  Wear  models,  however,  can  model 
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the  continuous  effect  of  deterioration  over  time.  Wear  models  were  introduced  histor¬ 
ically  later  than  shock  models  by  Esary  et  al.  [54] ,  which  proved  that  if  the  lifetime 
distribution  H(t)  of  a  device  is  subject  to  shocks  governed  by  a  Poisson  process  with 
the  probability  Pk  of  surviving  k  shocks,  and  if  the  discrete  distribution  is  increasing 
failure  rate  average  (IFRA),  then  II (t)  is  also  IFRA.  The  class  of  IFRA  life  distrib¬ 
utions  plays  a  fundamental  role  in  reliability  theory.  A  distribution  function  H  with 
survival  function  H  —  1  —  H  is  said  to  have  an  increasing  failure  rate  average  if 
H (0)  =  0  and  H~l{t )  is  decreasing  for  t  >  0.  Singpurwalla  [119]  provides  a  survey 
of  models  that  use  the  effects  of  a  common  environment  as  a  basis  for  generating 
dependent  lifetimes  for  components  in  a  system. 

Qinlar  [45]  introduced  a  Markov  additive  process  (MAP)  to  describe  the  failure 
mechanism  due  to  wear.  Markov-additive  processes  (MAPs)  are  a  class  of  Markov 
processes  which  have  important  applications.  A  MAP  {(A"(t),  J(t))  :  t  >  0}  is  a  bi¬ 
variate  Markov  process  whose  transition  probability  measure  is  translation  invariant 
in  the  additive  component  X(t).  Here  X{t )  is  any  independent  CTMC  from  which 
J (t )  is  an  additive  functional.  Qinlar  uses  the  unique  structure  of  the  MAP  to  prove 
that,  given  a  gamma  process  (a  stochastic  process  with  independent  increments)  with 
a  shape  parameter  that  is  a  function  of  Brownian  motion,  the  resulting  lifetime  is 
distributed  according  to  the  Weibull  distribution.  Kharoufeh  [66]  presents  a  compact 
transform  expression  for  the  failure  distribution  for  wear  processes  of  a  component 
degrading  according  to  a  Markovian  environment  inducing  state-dependent  continu¬ 
ous  linear  wear.  He  accomplishes  this  by  using  the  properties  of  a  MAP  and  assuming 
the  wear  process  to  be  temporally  homogeneous  and  that  the  environmental  process 
has  a  finite  state  space. 

Though  less  discussed  when  compared  to  individual  shock  and  wear  models, 
compound  damage  models  consider  the  combined  effect  of  wear  and  shocks  on  the 
lifetime  of  a  component.  Only  the  model  of  Kharoufeh  et  al.  [67]  considers  the 
deteriorative  effect  of  linear  wear  and  random  shocks  simultaneously.  Qinlar  [45] 
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derives  the  failure  time  distributions  and  associated  properties  for  models  subject  to 
continuous  wear  and  shocks  using  a  compound  Poisson  process.  However,  his  results 
are  complex  and  do  not  provide  a  computation  of  availability  measures.  Klutke  et 
al.  [71]  presents  availability  measures  for  a  component  with  random  inspection  times 
and  wear  rates.  Klutke  and  Yang  [70]  study  components  that  deteriorate  due  to  both 
shocks  and  graceful  degradation  with  periodic  inspections,  but  do  not  consider  the 
effect  of  wear  and  shocks  occurring  simultaneously.  Using  regenerative  arguments, 
they  derive  an  expression  for  the  limiting  average  availability.  Kiessler  et  al.  [68] 
studied  inspected  components  with  non-self- announcing  failures,  where  the  rate  of 
deterioration  is  governed  by  a  Markov  chain.  They  assumed  environmental  exposure 
changes  randomly  over  time  and  employed  a  Fourier  series  expansion  to  compute  the 
lifetime  distribution  and  availability  when  the  component  is  inspected  according  to 
a  periodic  inspection  policy. 

Taken  together,  the  shock,  wear  and  compound  damage  models  discussed 
herein  provide  a  number  of  mechanisms  useful  for  modeling  the  deterioration  of  a 
variety  of  systems.  However  degradation  models  provide  only  a  means  to  understand 
when  components  will  fail  and  do  not  provide  a  mechanism  to  maintain  them.  To 
that  end,  the  next  section  will  survey  optimal  maintenance  primarily  in  the  context 
of  optimal  replacement  models  which  maintain  a  system  through  replacing  failed 
components  according  to  a  wide  range  of  policies. 

2.3  Optimal  Replacement  Models 

Optimal  maintenance  models  have  received  a  great  deal  of  interest  in  the  last 
forty  years,  and  numerous  surveys  of  optimal  maintenance  have  been  contributed  to 
the  operations  research  literature.  The  first  was  completed  by  McCall  [87]  in  1965, 
which  covered  early  preventive  maintenance  models  and  classified  various  preventive 
maintenance  policies  by  their  common  characteristics.  That  same  year,  the  founda¬ 
tional  text  of  Barlow  and  Proschan  [29]  was  published.  Pierskalla  and  Voclker  [104] 
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and  Osaki  and  Nakagawa  [102]  provided  an  update  by  surveying  the  optimal  main¬ 
tenance  research  accomplished  in  the  11  years  elapsed  since  the  publication  of  [87]. 
These  papers  were  followed  by  the  survey  of  Sherif  and  Smith  [117]  in  1981  that  cat¬ 
egorized  524  different  optimal  maintenance  papers.  This  survey  was  followed  eight 
years  later  by  Valdez-Flores  and  Feldman  [131],  who  covered  optimal  maintenance 
research  since  [104].  More  recent  surveys  were  presented  by  Clio  and  Parlar  [43]  in 
1991,  Murdock  [93]  in  1995,  Dekker  [51]  in  1996  and  Wang  [138]  in  2002. 

Optimal  maintenance  models  can  be  classified  by  the  length  of  planning  hori¬ 
zon  for  the  problem,  the  degree  of  effectiveness  of  each  maintenance  action,  the 
optimization  criteria,  and  the  nature  of  the  degradation  process  and  state  space  that 
characterize  the  component  [131].  A  key  distinguishing  feature  among  these  models 
is  the  degree  to  which  a  component  is  repaired.  A  particular  advantage  of  perfect 
repair  models  is  that  they  are  able  to  take  advantage  of  renewal  theory,  as  presented 
in  [30]  and  [111].  This  section  will  briefly  survey  non-replacement  repair  models 
and  will  then  emphasize  and  survey  replacement  models  presented  in  the  literature. 
Replacement  or  renewal  refers  to  maintenance  models  in  which  a  given  component 
is  placed  into  a  state  after  repair  equal  to  the  condition  the  component  was  in  at  the 
start  of  the  process. 

The  first  alternative  to  replacement  or  perfect  repair  is  termed  minimal  repair. 
Minimal  repair  was  introduced  by  Barlow  and  Hunter  [28]  and  refers  to  models 
in  which  maintenance  actions  return  a  component  to  a  state  that  is  stochastically 
identical  to  the  state  just  prior  to  failure.  Thus  with  minimal  repair  the  failure  rate 
is  identical  before  and  after  a  repair  operation.  Minimal  repair  is  often  motivated 
by  considering  the  increased  cost  of  replacing  a  component,  when  compared  to  that 
of  a  simple  repair.  Several  minimal  repair  models  are  discussed  in  the  literature 
[113,  18,  118]. 

As  a  combination  of  perfect  and  minimal  repair,  imperfect  repair  is  incorpo¬ 
rated  into  models  in  which  maintenance  does  not  renew  a  component  to  its  original 
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state.  First  presented  by  Brown  and  Proschan  [36],  imperfect  repair  is  surveyed 
later  by  Pham  and  Wang  in  [139].  The  associated  mathematics  of  imperfect  repair 
is  more  involved,  as  standard  renewal  theory  does  not  apply.  Several  papers  present¬ 
ing  results  for  imperfect  repair  are  surveyed.  Recent  research  considering  imperfect 
repair  includes  Bruns  [37],  who  studies  a  repairable  component  with  Markovian  de¬ 
terioration  and  imperfect  repair  options  and  presents  optimal  strategies  to  minimize 
long-run  costs.  Biswas  et  al.  [33]  model  a  component  which  is  maintained  through 
a  periodic  inspections  with  maintenance  accomplished  by  a  combination  of  imper¬ 
fect  and  perfect  repair.  They  combine  the  maintenance  mechanisms  by  replacing  a 
failed  component  after  a  fixed  number  of  imperfect-repairs.  Cha  et  al.  [41]  compare 
steady-state  availabilities  of  two  different  components  subject  to  imperfect  repair 
policies.  They  make  comparisons  of  the  steady-state  availability  based  on  imperfect 
repair  policies.  Kijirna  and  Nakagawa  [69]  present  replacement  policies  for  a  shock 
model  with  imperfect  preventive  maintenance. 

Replacement  models  can  be  classified  as  either  preventive  or  corrective.  Cor¬ 
rective  maintenance  considers  all  maintenance  actions  performed  in  response  to  a 
component  failure.  Preventative  maintenance  is  concerned  with  actions  performed 
to  decrease  the  probability  of  unplanned  failure.  In  order  for  preventive  maintenance 
to  be  viable,  the  cost  of  a  preventive  repair  must  be  less  than  the  cost  of  repairing 
a  failure. 

There  are  two  primary  types  of  preventive  maintenance  policies:  block  replace¬ 
ment  and  age  replacement.  Under  a  block  replacement  policy,  maintenance  actions 
are  performed  at  fixed  and  deterministic  time  intervals.  Under  an  age  replacement 
preventive  maintenance  policy,  maintenance  is  initiated  when  a  component  has  accu¬ 
mulated  a  certain  amount  of  operational  time.  The  survey  papers  mentioned  before, 
especially  [87]  and  [104],  cover  many  variations  of  preventive  maintenance  models, 
with  an  emphasis  on  single  unit  components.  Research  in  the  area  of  age  replacement 
policies  is  primarily  concerned  with  the  optimal  operational  time  a  component  is  al- 
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lowed  to  accumulate  before  initiating  a  maintenance  action.  Common  assumptions 
to  accommodate  the  use  of  basic  renewal  theory  include  that  the  time  required  to 
replace  a  failed  component  is  considered  negligible  and  that  a  component  is  repaired 
perfectly  (i.e.  repaired  to  a  new  state).  Numerous  cost  models  have  been  presented 
to  determine  the  optimal  replacement  age,  most  notably  by  Ascher  and  Feingold  [19], 
A ven  and  Bergman  [26]  and  Nachlas  [94], 

Block  replacement  is  concerned  with  implementing  a  policy  of  maintenance  at 
regularly  spaced  intervals  without  regard  to  the  age  of  the  component.  Block  replace¬ 
ment  problems  are  popular,  partially  due  to  the  fact  that  they  are  often  easier  to 
implement  mathematically  and  operationally.  Barlow  and  Proschan  [30]  proved  that 
a  block  replacement  policy  results  in  less  failures  than  an  age  replacement  policy, 
even  though  block  replacement  requires  more  replacements  and  consequently  incurs 
higher  cost.  Wortman  et  al.  [144]  determined  that  deterministic  block  replacement  is 
preferable  to  random  replacement.  The  literature  discusses  many  preventive  mainte¬ 
nance  variations  of  block  replacement,  many  with  the  purpose  of  presenting  policies 
that  decrease  the  cost  of  frequent  replacement  or  simplify  the  mathematical  expres¬ 
sions  incurred  by  block  replacement.  For  example,  Barlow  and  Hunter  [28]  consider 
a  minimal-repair  preventive  maintenance  model  with  block  replacement.  Shaked 
and  Zhu  [114]  summarize  and  survey  various  preventive  maintenance  models  that 
use  renewal  theory  to  develop  the  mathematical  framework  for  block  replacement 
policies. 

Models  have  also  been  presented  that  combine  preventive  and  corrective  main¬ 
tenance,  such  as  the  work  by  Nachlas  and  Cassady  [95],  who  seek  to  determine  an 
optimal  balance  between  preventive  maintenance  and  corrective  maintenance.  Some 
papers  present  and  analyze  both  corrective  and  preventive  maintenance  as  potential 
maintenance  strategies,  while  placing  emphasis  on  a  particular  recommendation. 

While  block  replacement  models  suffer  from  the  potential  to  prematurely  re¬ 
place  components  and  age  replacement  policies  require  knowledge  of  the  cumulative 
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operational  time  for  the  component  under  consideration,  inspection  models  are  ap¬ 
plicable  when  the  state  of  a  component  is  known  only  during  an  inspection.  Many 
inspection  models  use  a  control  limit  policy,  which  determines  if  a  component  should 
be  replaced  upon  inspection  if  a  failure  or  cumulative  degradation  level  is  detected. 
Control  limit  policies  often  use  the  state  of  the  component  at  a  current  inspection 
to  determine  the  optimal  time  until  the  next  inspection.  The  potential  disadvantage 
of  inspection  policies  is  that  the  component  under  consideration  remains  in  a  failed 
state  until  the  next  inspection,  and  the  cost  of  a  component  in  a  down  state  must 
be  considered. 

One  of  the  first  inspection  models  presented  was  by  Luss  [83]  in  1976,  where 
the  state  of  a  component  is  modeled  by  several  levels  of  deterioration  and  the  cost  of 
replacement  after  a  failure  is  higher  than  before.  Zuckerman  [146]  reviews  inspection 
models  and  presents  an  optimal  inspection  interval  and  replacement  rule  with  consid¬ 
eration  for  the  associated  costs  when  degradation  is  due  to  a  Poisson  shock  process. 
The  research  presented  by  Wortman  et  al.  [144]  evaluates  two  different  inspection 
policies  for  non-self  announcing  failures  that  occur  as  a  result  of  a  Poisson  shock 
process.  Abdel-Hameed  [10]  studies  reliability  measures  when  degradation  occurs 
according  to  a  Markov  process  and  inspection  intervals  are  fixed.  More  recently, 
Vaurio  [136]  presents  an  inspection  model  with  replacements  occurring  according  to 
an  age  replacement  policy. 

Optimal  maintenance  models  can  be  classified  as  concerned  with  the  mainte¬ 
nance  of  either  a  single  unit  or  of  multiple  units.  Barlow  and  Proschan  [29]  introduce 
many  single-unit  models  and  authored  the  classic  textbook  used  for  the  mathemat¬ 
ical  theory  of  reliability.  Nakagawa  and  Osaki  [97]  followed  this  with  an  optimal 
replacement  model  that  considers  a  limited  number  of  spares.  Multiple-unit  main¬ 
tenance  models  often  suffer  from  a  state-space  explosion  problem.  While  much  has 
been  demonstrated  concerning  multiple  unit  maintenance  models  using  simulation, 
recent  research  uses  dynamic  programming  and  other  techniques  to  reduce  the  com- 
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plexity  incurred  by  an  expanding  state  size.  Ben  Ari  and  Gal  [31]  use  a  dynamic 
programming  approach  to  find  an  optimal  replacement  policy  for  a  multi-component 
system.  Flynn  et  al.  [56]  present  an  optimal  replacement  model  using  a  stochastic 
dynamic  program  with  the  goal  of  finding  the  optimal  balance  between  the  cost  of 
component  replacement  and  the  cost  of  component  failure.  Branch-and-bound  pro¬ 
vides  an  alternative  to  dynamic  programming  for  solving  multi-component  reliability 
problems,  as  presented  by  Chung  and  Flynn  [44], 

Availability  is  commonly  used  as  a  performance  measure  in  optimal  mainte¬ 
nance  models,  and  is  broadly  defined  as  the  fraction  of  time  a  component  is  in  an 
operational  state  over  the  total  time.  The  optimization  of  availability  is  both  an 
active  area  of  current  research  and  a  fertile  ground  for  future  work.  Maintaining 
a  component  for  optimum  availability  is  a  key  concern  in  many  military  and  op¬ 
erational  settings.  The  following  is  a  survey  of  the  optimal  maintenance  literature 
concerned  with  optimizing  and  measuring  availability  for  stochastically  degrading 
components. 

Wortman  and  Klutke  [145]  examine  the  availability  of  a  maintained  component, 
where  the  rate  of  deterioration  is  governed  by  an  random  environment.  With  the  ob¬ 
jective  of  exploring  the  influence  of  random  deterioration  on  component  availability, 
they  provide  a  result  that  exposes  the  relationship  between  the  remaining  lifetime, 
environment,  and  repairs.  Furthermore,  the  authors  develop  simple  bounds  that  can 
be  used  to  choose  inspection  rates  that  guarantee  a  specified  level  of  availability. 
The  principal  result  requires  no  specific  distributional  assumptions. 

Early  availability  models  include  Ahmed  and  Schenk  [15]  and  Srinivasan  and 
Ramachandra  [126].  In  1978  Ahmed  and  Schenk  [15]  used  optimal  control  theory  to 
obtain  the  optimal  policies  for  a  component  with  variable  failure  and  repair  rates, 
as  well  as  variable  maintenance  and  downtime  costs.  Their  optimization  model  con¬ 
sidered  the  cost  of  repair,  as  well  as  the  downtime  cost.  They  used  their  method  to 
consider  both  static  and  dynamic  optimal  maintenance  policies  and  note  that  the 
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optimal  maintenance  policy  is  highly  dependent  on  the  failure  rates  of  the  individual 
system  components.  Srinivasan  and  Ramachandra  [126]  applied  Kalman  filtering 
techniques  to  reliability  smoothing  and  used  these  to  compute  instantaneous  avail¬ 
ability  as  a  function  of  time.  From  this  they  were  able  to  approximate  the  sampling 
interval  and  the  measurement  variance  for  a  specific  maintenance  policy.  They  also 
proved  that  estimates  of  point  availability  converge  in  the  steady-state. 

Other  optimization  problems  have  been  presented  more  recently.  Rcineke  et 
al.  [109]  calculate  an  age  replacement  time  in  order  to  maximize  interval  availabil¬ 
ity.  They  extend  this  with  models  that  optimize  either  the  limiting  availability  or 
the  cost  rate  of  complex  components.  Cassady  et  al.  [39]  present  a  model  that 
approximates  the  replacement  age  in  order  to  maximize  the  average  availability  us¬ 
ing  experimental  design  and  regression  analysis.  More  recently  Cassady,  Pohl,  and 
Jin  [40]  presented  an  availability  optimization  model  that  applies  to  a  general  class 
of  two-state  repairable  components  by  determining  a  set  of  availability  importance 
measures.  Through  five  examples,  they  demonstrate  that  focusing  on  the  reduction 
of  component  failures  provides  greater  benefit  than  increasing  the  speed  of  equipment 
repair.  The  authors  then  established  a  set  of  three  optimization  models  that  cap¬ 
ture  the  trade-offs  between  improving  availability  performance  and  the  investments 
required  to  achieve  that  improvement. 

This  chapter  has  covered  a  broad  survey  of  the  literature  of  optimal  mainte¬ 
nance  relevant  to  this  thesis.  First,  the  history  of  optimal  maintenance  was  dis¬ 
cussed,  followed  by  an  overview  of  the  relevant  literature  of  three  key  degradation 
measures:  degradation  due  to  shocks,  degradation  due  to  wear,  and  degradation  due 
to  compound  models  combining  the  effects  of  shocks  and  wear.  Then  the  defining 
characteristics  of  optimal  maintenance  models  were  discussed  with  a  primary  em¬ 
phasis  on  replacement  models.  Availability  was  presented  as  a  useful  measure  to 
maximize  in  optimal  maintenance  models,  and  several  papers  were  reviewed  with 
associated  models  that  optimize  availability.  The  next  chapter  will  focus  on  the  par- 
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ticular  degradation  model  used  in  this  thesis:  a  compound  model  considering  wear 
due  to  a  random  environment  and  shocks  arriving  according  to  a  Poisson  process 
with  a  general  damage  distribution  and  an  associated  cost  function  that  considers 
costs  of  replacement,  downtime,  and  inspections. 
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3.  Mathematical  Model  and  Dynamics 

In  order  to  optimize  the  limiting  average  availability  by  selecting  an  appro¬ 
priate  inter-inspection  interval,  it  is  necessary  to  obtain  explicit  representations  for 
both  the  availability  and  cost  measures.  The  former  is  provided  by  Kharoufeh  et 
al.  [67] ,  who  derive  the  lifetime  distribution  as  well  as  the  limiting  average  availabil¬ 
ity  for  a  component  subject  to  linear  wear  and  random  shocks.  In  their  paper,  the 
authors  explicitly  derive  the  system’s  lifetime  distribution  and  mean  time  to  failure, 
as  well  as  an  equation  for  the  limiting  average  availability.  Due  to  the  relevance  of 
their  results  to  this  research,  their  main  results  are  briefly  reviewed. 

3.1  Stochastic  Degradation  Model 

Consider  a  component  evolving  in  a  random  environment  that  transitions  ac¬ 
cording  to  a  Markov  chain  in  continuous  time.  Degradation  is  due  to  the  combina¬ 
tion  of  two  damage  mechanisms:  environment-dependent  linear  wear  and  random 
shocks.  Degradation  accumulates  monotonically  until  an  inspection  detects  cumula¬ 
tive  damage  in  excess  of  a  fixed  threshold,  at  which  point  the  component  is  instantly 
replaced  with  one  in  new  condition.  The  maintenance  policy  is  to  inspect  the  cumu¬ 
lative  degradation  of  the  system  periodically  at  a  fixed  interval  of  length  r  and  to 
replace  the  component  if  degradation  is  found  to  be  beyond  a  fixed  threshold  value. 

Assume  the  environment  is  composed  of  a  finite  set  of  states,  each  correspond¬ 
ing  to  a  particular  constant  linear  wear  rate.  All  wear  rates  are  assumed  to  be 
positive,  which  implies  that  every  environmental  state  must  continuously  induce 
some  finite  linear  wear  per  unit  time  on  the  component.  The  second  degradation 
mechanism  considers  random  shocks  with  a  random  arrival  rate  and  a  single  general 
distribution  for  the  damage  caused  by  each  shock.  Shocks  are  assumed  to  arrive 
according  to  a  Poisson  process  with  rate  parameter,  A  where  A  is  a  positive  scalar. 
The  overall  effect  of  each  shock  is  considered  to  be  small  relative  to  the  degradation 
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process,  but  the  cumulative  effect  of  many  shocks  can  be  significant.  It  is  also  possi¬ 
ble  that  shocks,  unlike  wear,  may  cause  zero  damage.  The  wear  process  is  assumed 
to  be  stochastically  independent  of  the  shock  process. 

Component  inspections  are  assumed  to  be  instantaneous  and  perfect,  in  that 
there  is  no  measurement  error  in  determining  the  state  of  the  system.  Failures  are 
hidden  in  that  the  state  of  the  system  can  only  be  discerned  during  an  inspection. 
A  soft  failure  occurs  when  the  total  cumulative  degradation  exceeds  some  fixed 
threshold  value,  x.  The  time  at  which  this  threshold  is  crossed  can  be  considered 
the  component’s  useful  service  life.  Immediately  upon  inspection,  a  failed  component 
is  instantly  replaced  by  a  new  identical  component.  Partial  repairs  are  not  allowed. 

The  environment-dependent  rate  of  wear  is  modulated  by  a  continuous-time, 
stationary,  ergodic  continuous-time  Markov  chain  (CTMC),  Z,  where  {Zt  :  t  >  0}. 
Further,  the  finite  set  of  environmental  states  is  S  =  {1,2,...,/},  where  l  is  the 
dimension  of  the  state-space.  Each  state  has  a  corresponding  linear  wear  rate.  Thus 
Z  forms  a  finite,  continuous-time  Markov  chain  that  is  assumed  to  fully  characterize 
the  ambient  environment.  The  transitions  are  governed  by  an  infinitesimal  generator 
matrix  Q  with  an  initial  distribution  at.  The  transition  probability  functions  for  Z 
are  denoted  by 

Vij(t)  =  Pr  {Zt  =  j  |  Z0  =  i},  i,j  e  S. 

For  each  state  j  E  S  there  is  a  corresponding  positive  wear  rate  r(j).  Together, 
the  wear  rates  form  a  vector  r  of  dimension  /.  These  wear  rates  are  placed  along 
the  diagonal  of  a  matrix  of  zeros  to  form  the  diagonal  matrix  =  diag{r}.  The 
ratio  of  the  damage  threshold,  x,  to  the  minimum  wear  rate  is  denoted  by  A,  and 
is  considered  the  maximum  component  lifetime.  This  is  because  A  represents  the 
time  to  failure  if  the  component  were  subject  to  only  the  most  benign  damaging 
environment  until  failure  and  endured  no  damage  from  shocks.  A  is  formally  defined 
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by 


A  = 


x 


min{r(i)  :  i  G  S} 
The  total  accumulated  wear  Wt  at  time  t  is  defined  by 


(3.1) 


W,=  /  r(Zu)du. 

Jo 


Shocks  arrive  according  to  a  Poisson  process,  {Nt  :  t  >  0},  with  rate  A. 
The  sequence  of  independent  and  identically  distributed  (i.i.d.)  random  variables, 
{Yn  :  n  —  1,2, represents  the  individual  random  shock  magnitude  each  with 
cumulative  distribution  function  (c.d.f.)  FY.  The  Laplace-Sticltjes  transform  Fy(u) 
of  Fy(y)  is  given  by, 

poo 

Fy(u)=  /  e~uydFY(y), 

Jo 

and  the  diagonal  matrix  F D{u)  is  defined  to  be 

Fd{u)  =  diag{FY(u)}. 

The  total  accumulated  damage  j3t  at  time  t  is 

Nt 

n=0 

and  the  total  degradation  Xt  at  time  t,  can  then  be  expressed  as 

Xt  =  Wt  +  (Jt. 


The  time  required  for  Xt  to  first  reach  level  x  is  denoted  by  the  nonnegative 
random  variable  Tx\  i.e., 

Tx  =  inf {t  :  Xt  >  x}. 
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Since  no  partial  repairs  are  possible,  the  degradation  process  {Xt  :  t  >  0}  monotoni- 
cally  increases  until  the  first  inspection  after  Tx  or  the  start  of  the  next  replacement 
interval.  Furthermore,  the  bivariate  stochastic  process  {(X(t),Z(t))  :  t  >  0}  fully 
characterizes  the  state  of  the  system  at  time  t.  Since  Xt  is  monotonically  increasing 
and  all  wear  rates  are  positive,  the  event  in  which  total  cumulative  damage  is  below 
the  threshold  is  equivalent  to  the  event  in  which  total  time  elapsed  is  less  than  the 
time  of  failure,  or 

{Xt  <  x}  <=3-  {Tx  >  t}.  (3.2) 

The  process  describing  the  status  of  the  component  currently  in  use  is  the 
stochastic  process  (V’(t)  :  t  >  0}  given  by 

(  1  if  Xt  <  x 
^{t)  =  l 

[0  if  Xt  >  x 

Thus,  describes  the  state  of  the  system  at  time  t  and  is  recognized  as  a  right 
continuous  up/down  stochastic  process.  The  nth  component  lifetime  is  denoted  by 
Ln  with  a  corresponding  time  of  failure  designated  by  Fn.  Since  the  state  of  the 
system  is  only  discerned  at  each  inspection,  the  nth  replacement  epoch  Rn  occurs  at 
the  first  inspection  after  the  nth  failure  time,  or 

Rn  =  min{h  :  kr  >  Fn},  n  —  0, 1,  2, ... . 
with  Rq  dehned  to  be  0. 

The  entire  process  is  graphically  depicted  through  the  sample  paths  displayed 
in  Figure  3.1.  Without  loss  of  generality,  assume  that  the  environment  considered 
in  Figure  3.1  can  be  adequately  described  by  three  states.  In  this  example,  the 
first  environment  has  the  most  corrosive  associated  linear  wear  rate,  while  the  third 
environmental  state  has  the  most  benign  rate  of  wear.  Besides  demonstrating  the 
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Figure  3.1  Sample  paths  of  Xt  and  Zt. 

cumulative  deterioration  process  as  a  combination  of  linear  wear  and  random  shocks 
subject  to  a  control  limit  policy,  the  sample  path  analysis  demonstrates  that  succes¬ 
sive  lifetimes  do  not  form  an  i.i.d.  sequence  of  random  variables  unless  they  begin 
their  lifetimes  in  identical  initial  environmental  states.  Moreover,  the  dependency  of 
failure  and  replacement  epochs  on  component  lifetimes  is  clear. 

Also  of  interest  is  the  embedded  process  representing  the  state  of  the  environ¬ 
ment  at  the  nth  replacement  epoch,  {£n  :  n  >  0}.  Since  there  are  a  finite  number  of 
states,  £n  forms  an  irreducible  embedded  discrete  time  Markov  chain  (DTMC)  on  S 
with  transition  probability  matrix  P  and  stationary  distribution  p.  An  important 
result  reported  in  [67]  is  that  {(£n,i?n)  :  n  >  0}  forms  a  Markov  renewal  process. 
This  is  formally  stated  in  Theorem  3.1  and  a  formal  proof  is  available  in  [67].  This 
theorem  is  necessary  in  order  to  use  an  availability  equation  that  requires  a  renewal 
process  since  the  sequences  {Fn  :  n  >  0}  and  {Rn  :  n  >  0}  do  not  form  a  set  of  i.i.d. 
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component  lifetimes  unless  the  replacement  epoch  occurs  when  the  environment  is 
in  the  same  initial  state. 


Theorem  3.1  The  bivariate  process  {(£n,  Rn)  :  n  >  0}  is  a  Markov  renewal  process, 
and  the  continuous-time  process  {^{t)  :  t  >  0}  is  Markov-regenerative  w.r.t  {(£n,  Rn)  : 
n  >  0} 


Let  E  [•]  denote  the  expectation  operator,  and  E.,  [•]  denote  the  conditional 
expectation  E  [•  |  Zq  —  i] .  The  limiting  average  availability  is  defined  by  Leernis  in 
[76]  as, 


A  =  lim  t  1  /  E  [' ip(w)\  dw.  (3.3) 

t^°o  JQ 

Theorem  3.2,  which  was  proved  by  Kiessler  et  al.  [68:704],  expresses  the  limiting 
average  availability  as  a  direct  consequence  of  the  renewal  reward  theorem.  The¬ 
orems  3.1  and  3.2  allow  for  the  limiting  average  availability  to  be  represented  as 
the  expected  time  to  the  first  failure  divided  by  the  expected  time  to  the  first  re¬ 
placement,  only  if  these  expected  values  are  conditioned  on  the  initial  state  of  the 
environment. 


Theorem  3.2  The  limiting  average  availability  is  given  by 

A(T )  =  lim  r1  f  EIV-W]  dw  =  (3.4) 

^°°  Jo  Ei= i  Pi [-Ri] 

where  pi  =  lim,woo  Pr  {^n  —  i}  with  i  —  1,2, ...  ,1. 

In  order  to  compute  equation  (3.4)  Kharoufch  et  al.  [67]  derived  explicit  representa¬ 
tions  for  Ej  [F\\  and  Ej  [i?i]  in  terms  of  Q,  A,  Fy(w),  R^,  x  and  r. 

The  distribution  function  of  the  total  random  lifetime  of  the  system  is  defined 
as 

G(x,t)  =  Pr  {Tx  <  t} 
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and  the  component  lifetime  distribution  conditioned  on  the  initial  state  i,  is  denoted 
as  Gi(-)  =  G(-  |  Zq  —  i).  The  conditional  distribution  function  can  be  written  as, 

Gi(x,  t )  =  Pr  {Tx  <  t  |  Zq  —  i}  ,  i  e  S. 

Kharoufeh  et  al.  [67]  show  that  the  failure  time  distribution  satisfies  a  system 
of  linear  first-order  partial  differential  equations.  They  solve  this  system  using  trans¬ 
form  methods,  yielding  the  Laplace-Stieltjes  transform  of  the  lifetime  distribution. 
For  further  details,  the  reader  is  referred  to  [67].  This  distribution  can  be  used  to 
determine  the  Laplace-Stieltjes  transform  of  the  conditional  and  unconditional  mean 
lifetime  with  respect  to  x  and  also  determine  the  conditional  expectation  of  the  first 
replacement  epoch. 

Component  conditional  and  unconditional  lifetime  distributions  are  stated  in 
Theorem  3.3.  In  this  theorem,  a.  denotes  the  initial  state  probability  vector,  e 
denotes  a  column  vector  of  ones,  and  e,  denotes  the  ith  column  in  an  identity  matrix 
of  dimension  /. 

Theorem  3.3  The  Laplace-Stieltjes  transform  of  the  unconditional  and  conditional 
first  lifetime  distributions,  G(x,t )  and  Gi(x,t),  with  respect  to  x  are,  respectively, 

G(u,  t)  —  1  —  a  exp  +  A(F d(u)  -  I)  -  «  Rp)  j  i j  e  (3.5) 

and 

Gi(u,t )  —  1  —  ej  exp  (^Q  +  A(F D{u)  -  I)  -  wRd))  e.  (3.6) 

Theorem  3.1  requires  that  the  expected  times  to  first  failure  and  first  replace¬ 
ment  must  be  conditioned  on  the  initial  environmental  state,  necessitating  an  ex¬ 
pression  for  the  transition  probability  matrix  P  for  the  DTMC  {£n  :  n  >  0}.  The 
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(i,  k)th  element  of  P  is 


Pi,k  =  ^2  7r*.fc(nr)  Pr  iRi  =  nT  \  zo  =  0  (3.7) 

n=l 

where,  7  is  the  number  of  inspection  intervals  required  to  exceed  the  maximum 
component  lifetime,  or 

7  =  min{n  >  1  :  nr  >  A}.  (3.8) 

Since  7  determines  the  number  of  terms  which  must  be  summed  according  to  equa¬ 
tion  (3.7),  its  size  is  very  important  to  the  overall  computation  speed.  For  use  in  a 
computational  setting,  the  ceiling  function  can  equivalently  be  used  to  define  7  as 

7  =  [A/r]  .  (3.9) 

The  probability  Pr{i?!  =  nr  |  Z0  =  i}  is  given  by 

Pr  {i?x  =  nr  \  Z0  =  i}  =  Gi(x,  nr)  —  G^x,  {n  —  l)r)  =  A i(x,  nr).  (3.10) 

Due  to  the  structure  of  the  process,  the  expected  time  of  the  first  failure  is 

Ej  [Fi]  =  E  [Tx  j  Z0  =  i] .  (3-11) 

By  applying  the  definition  of  expectation  and  employing  the  properties  of  Laplace 
transforms  to  (3.6),  the  Laplace-Stieltjes  transform  of  (3.11)  becomes 

E  [Tu\  =  or  (uRd  -  Q  -  A  (fd(u)  -  ijj  e  (3.12) 

with  the  conditional  representation, 

E.t  [Tu]  =  ej  (uRd  -  Q  -  A  (f  D(u)  -  ijj  e,  (3.13) 
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as  shown  in  [67]. 

The  sole  expression  remaining  to  determine  the  limiting  average  availability  is 
the  conditional  expected  time  for  the  first  replacement  to  occur.  By  proving  that  Tx 
is  bounded  above  by  x/ri,  where  ry  =  min{r(i)  :  %  G  S}  is  the  minimum  linear  wear 
rate,  and  by  conditioning,  Kharoufeh  et  al.  [67]  compute  E.,  [i^]  as, 


Ei  [i?i]  =  r  f  7  -  y ^Gj(x,  nr] 


(3.14) 


With  expressions  for  the  expected  time  of  the  first  failure,  and  the  expected 
time  for  the  first  replacement,  the  limiting  average  availability  can  now  be  computed. 
From  Theorem  3.2  and  the  various  components  defined  above,  the  limiting  average 
availability  can  be  expressed  as 


A(t)  = 

Eli  Pi®i[Ri] 

El;=i Pi  | -  Q  -  A  (fb(m)  -  i))  e|^) 
Y!i=iPi'r  (7  -  EEo  {«_1Gj(M,nr)}) 


(3.15) 


where  C  1  {•}  is  the  inverse  Laplace  transform  operator. 


3.2  Cost  Function  Development 

The  objective  of  this  research  is  to  maximize  the  limiting  average  availability  as 
represented  in  equation  (3.15)  by  selection  of  the  inter-inspection  interval  (r)  while 
remaining  within  an  arbitrary  budget  constraint  B  and  incorporating  a  realistic  cost 
scheme.  Further,  for  this  computation  to  be  operationally  useful,  it  is  necessary  to 
compute  the  maximum  limiting  average  availability  with  minimal  computation  time. 
This  optimization  will  be  accomplished  by  defining  a  formal  nonlinear  program  us¬ 
ing  equation  (3.15)  as  the  objective  function,  constrained  by  cost  and  boundaries 
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imposed  by  problem  structure.  Computational  performance  improvements  are  dis¬ 
cussed  in  chapters  4  and  5. 

In  order  to  maximize  availability  for  a  particular  inspection  interval,  it  is  nec¬ 
essary  to  consider  a  set  of  costs  relevant  to  industrial  or  military  applications  that 
remain  within  an  associated  budget,  B.  The  renewal  reward  theorem  [111]  allows  for 
consideration  of  the  long-run  costs.  Three  long-run  costs  per  each  cycle  are  consid¬ 
ered:  (1)  fixed  cost  of  replacing  a  failed  component  Cr,  (2)  cost  of  downtime  Cr,  and 
(3)  per-inspection  cost  Cj.  The  expected  cycle  time  is  equivalent  to  the  expected 
time  to  the  first  replacement  or, 

i 

E  Pi  E*  iRl\  ■ 

i= 0 

When  divided  by  the  expected  cycle  time,  the  sum  of  these  three  long-run  costs  forms 
the  long-run  expected  cost  per  cycle,  which  is  limited  by  the  availability  improvement 
budget  constraint,  B.  The  long-run  cost  per  cycle  is, 

CR  +  C!  +  Cd 
E  [R^ 

The  replacement  cost,  Cr,  is  understood  as  the  cost  to  replace  a  failed  com¬ 
ponent.  Since  a  replacement  marks  the  end  of  a  cycle,  this  cost  occurs  exactly  once 
per  cycle  and  is  therefore  a  fixed  quantity,  not  dependent  on  r.  In  most  industrial 
settings  this  cost  is  much  higher  than  the  per-inspection  cost  and  often  higher  than 
the  cost  per  unit  of  downtime.  Furthermore,  this  cost  represents  an  aggregate  of  the 
labor  and  material  costs  involved  in  replacing  a  component. 

A  cost  vital  to  industry  and  defense  is  the  cost  incurred  per  unit  time  for  com¬ 
ponent  to  be  in  a  non-operational  state.  In  an  industrial  setting,  this  cost  per  unit 
of  downtime  could  be  the  revenue  lost  from  a  machine  producing  needed  inventory. 
In  a  military  setting,  this  downtime  could  be  the  operational  cost  of  not  having  the 
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use  of  a  mission-critical  protection  mechanism.  The  cost  per  unit  of  downtime,  cd 
varies  widely  depending  on  the  criticality  of  the  component.  The  total  long-run  cost 
of  downtime  Cd  in  a  particular  cycle  is  the  cost  per  unit  of  downtime  multiplied  by 
the  expected  elapsed  time  between  a  system  failure  and  the  start  of  the  next  cycle. 
Thus,  Cd  is  given, 


Cd  —  cd  E  [Ri  —  F\  —  cd  E  [.Ri]  —  cd  E  [T\] . 


Therefore,  the  greater  the  expected  elapsed  time  between  a  hidden  failure  and  the 
next  inspection,  the  greater  the  long-run  cost  of  downtime.  A  smaller  r  presumably 
ensures  a  smaller  downtime  cost  as  more  frequent  inspections  decrease  the  expected 
time  from  a  failure  to  the  next  inspection. 

The  final  long-run  cost  considered  is  the  long-run  expected  cost  per  inspection 
C f,  which  occurs  only  when  the  cycle  time  is  large  enough,  or  if  Ri  >  r.  The  long- 
run  number  of  inspections  in  a  cycle  is  the  expected  length  of  the  cycle  divided  by 
the  deterministic  inter-inspection  interval  length.  Since  partial  inspections  have  no 
meaning,  the  remainder  of  this  quotient  must  be  excluded.  If  a  single  inspection  costs 
cj,  then  the  expected  total  cost  due  to  inspections  in  a  given  cycle  is  cj  multiplied 
by  the  expected  number  of  inspections  within  a  given  cycle.  Mathematically,  the 
inspection  cost  can  be  expressed  using  the  floor  function  |_-J  as 


C!  =  Cl 


E  [Ri] 


The  longer  the  periodic  inspection  interval,  the  lower  the  long-run  cost  clue  to  inspec¬ 
tions.  Conversely,  frequent  inspections  cause  high  inspection  costs.  In  this  manner, 
the  per-inspection  cost  penalizes  frequent  inspections.  In  an  industrial  setting,  a 
single  per-inspection  cost  is  often  far  lower  than  the  cost  of  down-time  and  the  cost 
of  replacement. 
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Since  successive  unit  lifetimes  form  an  i.i.d.  sequence  of  random  variables  only 
if  they  are  placed  into  service  in  identical  environmental  states,  the  conditional  ex¬ 
pectations  for  the  first  failure  and  first  replacement  must  be  multiplied  by  the  sta¬ 
tionary  probabilities  p,  and  summed  to  produce  the  unconditional  expected  times 
to  first  failure  and  first  replacement.  By  the  renewal  reward  theorem,  the  long-run 
cost  is  the  sum  each  of  the  expected  costs  per  cycle  divided  by  the  expected  cycle 
length,  which  is, 


Ei= 1  Pi  {Cr  +  Cr>  Ej  [i?i]  —  Cr  Ej  [Fi]  +  a 


ELiP.E.  [ft] 


(3.16) 


It  is  important  to  note  that  this  cost  function  is  not  convex,  since  the  inspection  cost 
is  not  a  continuous  function  with  respect  to  r.  Moreover,  since  equations  for  Ej  [iq] 
and  Ej  [f?i]  are  represented  only  in  the  form  of  Laplace  transforms,  equation  (3.16) 
is  not  guaranteed  to  be  smooth  or  continuous.  As  will  be  discussed  in  chapter  4,  this 
precludes  the  application  of  traditional  nonlinear  optimization  techniques. 


3.3  Mathematical  Programming  Formulation 

In  order  to  maximize  the  limiting  average  availability  as  represented  in  equa¬ 
tion  (3.15)  by  selection  of  the  inter-inspection  interval,  the  following  nonlinear  pro¬ 
gram  was  constructed.  The  objective  function  requires  x,  Rd,  Fo(m),  A  and  Q  to  be 
defined  a  priori  and  evaluates  the  limiting  average  availability  as  a  function  of  the 
decision  variable,  r,  the  deterministic  inspection  interval.  The  full  objective  function 
is, 

iPi  jw_1ef  (uRd  -  Q  -  A  (f d(u)  -  i))  ej 
ELi  Pi  T( 7  -  En=o£~‘  {w^GjKnr)  j) 
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which  can  be  written  more  compactly  as 


A(t) 


EUftE*  [-El] 


3. 3. 1  Constraints 


In  order  to  form  a  complete  optimization  problem,  bounds  on  r  were  computed 
for  the  limiting  average  availability  as  presented  by  the  stochastic  degradation  model. 
It  was  proved  by  [55]  that  Tx  is  bounded  by  A.  Moreover,  when  r  reaches  A,  the  value 
of  7  is  unity,  and  the  expected  time  for  the  first  replacement  epoch,  as  presented 
in  equation  (3.14),  is  also  A.  As  the  value  of  r  increases  past  A  by  some  positive 
valued  quantity  5,  the  value  of  7  remains  at  one  since  one  inspection  interval  is  large 
enough  to  cover  A  or, 


A 

A  +  h 


1,  V  <5  >  0. 


Therefore,  for  some  77  >  A  where  77  =  A  +  6,  equation  (3.14)  becomes 


E,  [Ri]  =  77 


7 


7-1 

y^Gi(x,nr) 


n= 0 


A  +  5. 


(3.17) 


Therefore,  the  denominator  of  A(t)  increases  proportionally  to  5  and  A(t)  ap¬ 
proaches  zero  since  the  numerator  is  constant  with  respect  to  r.  However,  inspecting 
the  system  after  a  failure  has  occurred  with  probability  one,  does  not  make  sense 
in  an  operational  context.  Thus,  the  upper  bound  of  r  is  set  to  A  and  E,  [Ah]  is 
consequently  less  than  or  equal  to  A. 

In  order  to  maximize  the  limiting  average  availability,  it  is  necessary  to  un¬ 
derstand  availability  measures  as  r  approaches  zero.  In  words,  inspecting  infinitely 
often  generates  an  availability  of  one.  This  is  understood  in  that,  as  r  gets  smaller, 
the  expected  time  until  the  next  replacement  approaches  the  expected  time  to  the 
next  failure  and  failures  cease  to  be  hidden  in  that  they  are  instantly  detected.  How¬ 
ever,  while  r  can  get  very  small,  it  will  always  be  positive  with  a  finite  B.  This  is 
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due  to  the  inversely  proportional  relationship  between  Cj  and  r.  Inspection  cost  ap¬ 
proaches  infinity  as  r  gets  small.  That  the  decision  variable  remains  finite  over  the 
interval  of  interest  is  an  important  property  for  several  optimization  methodologies. 

Even  though  both  the  objective  function  and  the  cost  constraint  must  be 
treated  as  black-box  functions,  both  satisfy  the  Lipschitz  condition,  which  is  neces¬ 
sary  to  prove  convergence  for  a  number  of  derivative-free  optimization  methods.  A 
function  /  :  /  — >  M  is  said  to  be  Lipschitz  continuous  if  there  exists  a  constant  K  >  0 
such  that,  |  f(x)  —  f(y) \  <  K\x  —  y\  for  all  x,  y  in  the  interval  /.  By  definition,  A{r) 
produces  values  in  the  region  [0,1].  As  established,  the  interval  of  M  for  this  problem 
is  restricted  to  (0,  A] .  Thus  for  any  T\  and  r2  in  (0,A], 

|A(ti)  -  A(t2)  |  <  1 

and  the  denominator  is  bounded  by, 


ri  -  r2|  <  A. 


Therefore,  with  the  Lipschitz  constant  0  <  K  <  oo, 

Tl, 72,(0,  a], 

Ti  _  r2  A 


the  objective  function  is  Lipschitz  continuous. 

By  the  constraints  imposed  on  r,  Ej  [i?x]  <  A.  For  a  finite  r,  Rn  >  Fn  and  the 
expected  cycle  time  will  be  positive  as  long  as  r  is  nonzero.  Therefore,  E;  [iih]  — Ej  [T1!] 
will  always  be  some  finite  number.  Moreover,  the  finite  budget  constraint  ensures  the 
long-run  cost  is  always  finite.  Thus,  the  cost  constraint  is  also  Lipschitz  continuous 
since 


|C'(ti)-C'(t2)| 


Ti  —  r2 


<K  Ti,t2g(0,A]. 
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3.3.2  Nonlinear  Programming  Formulation 

With  the  constraints  and  the  decision  variable  defined,  the  full  nonlinear  pro¬ 
gram  can  be  stated  as 

Maximize  A(t)  =  ^I:=1  Pi  E*  ^  (3.18) 

EL  Pi^ilRi] 

subject  to 

Ell  Pi  (cR  +  CD(  E,  [R,]  -  Ej  [F,])  + 

- - - t - -  <  B 

Ei=i  P‘  Ei  [-Ri] 

0  <  t  <  A 

This  chapter  has  summarized  the  results  of  the  stochastic  degradation  model 
presented  by  Kharoufeh  et  al.  [67],  in  which  degradation  was  the  result  of  random 
shocks  and  environment-dependent  linear  wear.  From  this  paper,  the  availability  and 
reliability  measures  were  presented  and  discussed.  Then,  a  cost  model  to  consider 
the  long-run  costs  using  the  renewal  reward  theorem  was  developed.  Constraints  on 
r  were  discussed  and  developed.  The  availability  and  reliability  measures  from  the 
stochastic  degradation  model  and  the  the  cost  constraint  function  were  combined  to 
formulate  a  nonlinear  program. 

However,  the  equations  for  Ej  [Fj]  and  E,  [i?i]  can  be  only  represented  us¬ 
ing  Laplace  transforms  and  must  consequently  be  numerically  inverted  to  compute. 
Moreover,  the  transition  probability  matrix  P  requires  computing  the  component 
lifetime  distribution  which  also  requires  the  numerical  Laplace  transform  inversion. 
This  presents  considerable  challenges  and  precludes  the  use  of  many  standard  op¬ 
timization  techniques.  The  next  chapter  will  discuss  a  solution  methodology  and 
optimization  strategies  to  solve  the  nonlinear  program  of  section  3.3. 
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4.  Solution  Procedures 


This  chapter  discusses  the  approach  used  to  approximately  solve  the  nonlinear 
programming  problem  presented  in  section  3.3.  First,  standard  solution  procedures 
for  typical  optimization  problems  are  discussed.  While  modifications  to  these  pro¬ 
cedures  for  nonsmooth  objective  functions  are  promising,  they  are  generally  not 
well-suited  for  an  objective  function  that  exists  only  in  the  form  of  a  Laplace  trans¬ 
form.  However,  solutions  may  be  obtained  by  using  generalized  pattern  search,  an 
approach  that  requires  neither  smooth  functions  nor  derivative  information  and  has 
proven  convergence  properties.  After  a  discussion  of  pattern  search,  an  algorithm 
directly  implementing  the  analytical  model  of  chapter  3  to  compute  the  limiting 
average  availability  A  is  presented,  followed  by  several  strategies  to  improve  compu¬ 
tational  performance.  These  result  in  an  improved  iterative  algorithm  for  maximizing 
A(t)  with  a  truncation  method  that  dramatically  decreases  the  computational  re¬ 
quirements  necessary  to  solve  the  original  optimization  problem.  Finally,  in  order 
to  better  improve  computational  performance,  two  of  the  most  often  computed  and 
consequently  most  computationally  expensive  algorithms  are  discussed:  the  matrix 
exponential  and  the  Laplace  transform  inversion  algorithm. 

The  optimization  problem  of  Section  3.3  seeks  to  compute  the  maximum  of  the 
limiting  average  availability  A(r)  on  the  interval  (0,  A]  by  selection  of  the  value  of  r 
that  maximizes  A(r),  denoted  by  t*.  Since  the  objective  and  cost  functions  are  not 
necessarily  smooth  or  even  continuous,  traditional  necessary  and  sufficient  conditions 
for  optimality  may  not  apply.  Since  the  problem  is  nonconvex  and  treated  as  a  black 
box ,  convergence  to  a  global  optimizer  cannot  be  guaranteed  [127].  Therefore,  the 
variable  r*  is  assumed  be  only  a  local  maximizer,  which  means  that,  for  some  £  >  0, 
A(t*)  >  A{t)  for  all  r  e  (0,  A]  satisfying  ||r  —  r*||  <  e. 

As  discussed  in  chapter  3,  the  objective  function  in  equation  (3.18)  involves 
a  Laplace  transform  and  contains  a  generally  defined  c.d.f.  This  adds  considerable 
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complexity  to  the  overall  optimization  problem,  since  numerical  methods  must  be 
used  to  compute  inverse  Laplace  transforms  and  derivative  information  is  not  avail¬ 
able.  Moreover,  the  objective  function  may  not  be  smooth,  since  no  such  assumption 
is  made  on  the  c.cl.f. 

Clearly,  an  analytical  solution  to  (3.18)  is  preferred,  in  which  elementary  cal¬ 
culus  yields  a  simple  expression  for  the  optimal  r.  It  is  sometimes  possible  to  apply 
partial  fraction  decomposition  to  a  Laplace  transform-based  function  in  order  to  fa¬ 
cilitate  analytic  inversion  using  Laplace  transform  tables.  However,  in  this  research, 
an  analytical  representation  for  neither  the  first  replacement  epoch  (E,  [Ri])  nor  the 
first  failure  epoch  (Ej  [E)])  was  attainable;  consequently,  numerical  methods  were 
required  to  compute  both  the  objective  function  and  the  nonlinear  cost  constraint. 

Many  numerical  algorithms  exist  for  solving  unconstrained  optimization  prob¬ 
lems  with  differentiable  functions,  including  gradient  descent,  Newton’s  method, 
quasi-Newton  methods,  and  conjugate  gradient  methods.  These  techniques  require 
an  initial  feasible  point  and  use  derivative  information  to  generate  a  sequence  of 
iterates  that  converge  to  a  point  satisfying  certain  optimality  conditions.  Extensions 
to  these  methods  for  constrained  problems  also  appear  throughout  the  literature 
(see  [32],  for  example).  If  derivatives  are  not  available,  a  common  practice  is  to  ap¬ 
proximate  them  with  finite  differences.  However,  their  use  still  assumes  that  deriv¬ 
atives  exist,  even  though  they  are  not  available.  Modifications  of  derivative-based 
methods  for  nonsmooth  functions  also  exist  [106,  103,  64,  63],  but  they  generally  re¬ 
quire  an  explicit  objective  function  to  exist  at  every  point  on  the  interval  of  concern 
in  the  real  domain. 

Heuristics,  such  as  local  search,  greedy  algorithms,  simulated  annealing,  ant 
colony  optimization,  tabu  search  and  genetic  algorithms  (e.g.,  see  [90]),  are  useful 
for  problems  in  which  other  methods  cannot  be  applied.  They  treat  the  objective 
function  as  a  black  box  and  are  useful  in  practice  for  improving  the  objective  function 
value;  however,  useful  convergence  properties  are  extremely  rare  [91],  and  the  random 
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nature  of  some  of  these  methods  would  undoubtedly  generate  trial  points  r  close  to 
zero,  incurring  unnecessary  computational  expense  due  to  the  associated  large  values 
of  7. 

Among  the  multitude  of  direct  search  methods,  three  other  approaches  were 
considered:  Nclder-Mead,  golden  section  search,  and  generalized  pattern  search 
(GPS).  Nelder-Mead  is  a  flexible  numerical  method  for  optimizing  multidimensional 
unconstrained  problems  and  belongs  to  the  more  general  class  of  direct  search  algo¬ 
rithms  [98],  Further,  Nelder-Mead  is  a  commonly  used  for  nonlinear  optimization 
problems  and  known  to  perform  well  when  the  objective  function  is  difficult  or  expen¬ 
sive  to  compute,  not  smooth,  or  when  function  evaluation  is  noisy  [101],  Nelder-Mead 
can  work  effectively  in  these  situations  because  it  only  uses  function  evaluations  to 
reach  a  solution.  Nelder-Mead  has  been  shown  to  converge  to  the  optimal  solution 
when  applied  to  strictly  convex  functions  in  one  and  two  dimensions  [74],  but  not  in 
the  general  case.  In  fact,  counter-examples  exist  [88]  that  show  that  Nelder-Mead 
cannot  even  guarantee  convergence  to  a  first-order  stationary  point. 

Golden  section  search  is  a  method  specifically  for  locating  the  minimum  point 
of  a  continuous  function  on  a  fixed  interval.  It  requires  no  information  about  the 
derivative  of  the  function.  Golden  section  search  starts  with  a  function  evaluation 
f(x)  at  some  point  x  in  the  larger  of  the  two  intervals  (a,  b )  or  (6,  c).  If  f(x)  <  f(b), 
then  x  replaces  the  midpoint  b,  and  b  becomes  an  end  point.  If  f(b)  <  f{x ),  then 
b  remains  the  midpoint  with  x  replacing  one  of  the  end  points.  Regardless  of  the 
outcome,  the  width  of  the  bracketing  interval  will  decrease.  The  procedure  is  re¬ 
peated  until  the  final  width  achieves  a  desired  tolerance.  If  the  new  test  point  x  is 
chosen  to  be  the  specific  proportion  of  (3  —  a/5)/2  (known  as  the  golden  section) 
along  the  larger  sub-interval,  measured  from  the  mid-point  b ,  then  function  values 
can  be  re-used,  and  the  number  of  function  values  to  achieve  a  desired  level  of  accu¬ 
racy  is  minimized.  However,  the  method  converges  slowly,  cannot  handle  nonlinear 
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constraints,  and  is  not  extendable  to  multiple  dimensions.  For  further  details,  the 
reader  is  referred  to  [105]. 

GPS  is  a  derivative-free  method  originally  introduced  by  Torczon  [130]  for 
unconstrained  problems,  who  proved  convergence  of  a  subsequence  of  iterates  to 
first-order  stationary  point.  It  has  known  convergence  properties  for  a  variety  of 
problem  classes,  even  when  the  objective  function  is  nonsmooth  [22],  Generalized 
pattern  search  methods  iteratively  search  a  set  of  points  around  the  current  iterative 
point  for  one  that  improves  the  objective  function  value.  Because  pattern  search 
methods  do  not  require  derivative  information  to  search  this  set  of  points,  they 
are  commonly  used  when  evaluation  of  the  objective  function  is  expensive  or  when 
accurate  approximation  of  derivatives  is  problematic.  Consequently,  pattern  search 
can  be  applied  to  the  problem  in  (3.18),  even  though  the  objective  function  requires 
inversions  of  several  Laplace  transforms. 

Due  to  the  shortcomings  of  the  other  methods  described  in  this  section,  GPS 
was  chosen  as  the  optimization  method  for  this  thesis.  The  next  section  describes 
the  algorithm  in  greater  detail,  and  in  a  manner  similar  to  that  of  [12]  and  [22], 
Convergence  results  are  discussed  shortly  thereafter. 

4-1  Generalized  Pattern  Search 

Consider  the  general  nonlinear  minimization  problem, 

min  f(x),  (4.1) 

x£il 

where  D  =  {x  e  M"  :  £  <  Arc  <  it},  /  :  Mn  — »  M  and  A  e  Qmxn  is  a  rational  matrix. 
Moreover,  £  and  u  are  the  lower  and  upper  bounds  of  the  constraints  where  £  and 
u  are  e  (Mm  n  {±oo}}  and  £  <  u. 

GPS  algorithms  generate  a  sequence  of  iterates  {re*,}  in  Mn  with  nonincreasing 
objective  function  values.  Each  iteration  is  divided  into  an  optional  SEARCH  step  and 
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a  local  poll  step.  Both  the  search  and  poll  steps  evaluate  points  on  a  mesh  in 
order  to  find  a  mesh  point  with  an  improved  function  value.  The  mesh  is  constructed 
as  a  lattice  of  points  in  Mn,  based  on  a  finite  set  of  directions  D  that  form  a  positive 
spanning  set  and  a  mesh  size  parameter  Ak  >  0  that  controls  the  fineness  of  the 
mesh.  In  this  case,  a  positive  spanning  set  refers  to  a  set  of  vectors  such  that  any 
vector  in  the  space  to  be  represented  by  a  nonnegative  linear  combination  of  the 
vectors  in  the  set. 

By  definition,  nonnegative  linear  combinations  of  the  elements  of  the  set  D 
span  Mn.  The  directions  that  form  D  can  be  arbitrarily  chosen  provided  that  each 
direction  dj  E  D,  j  =  1,  2, . . . ,  |T>|,  dj  =  G Zj,  where  G  G  Mnxn  is  a  nonsingular 
matrix  and  z3  E  Zn  is  an  integer  vector.  At  iteration  k ,  the  mesh  is  centered  around 
the  current  iterate  Xk  G  M"  and  its  fineness  is  parameterized  through  the  mesh  size 
parameter  A The  mesh  can  then  be  represented  as 

Mk  =  +  AfcD  z  :  z  E  zj1}  (4.2) 

where  Z+  is  the  set  of  nonnegative  integers.  Note  that  in  (4.2),  the  columns  of  the 
matrix  D  form  the  set  D. 

In  the  SEARCH  step,  GPS  can  evaluate  any  finite  set  of  mesh  points,  and  a  num¬ 
ber  of  strategies  exist  for  generating  trial  points,  including  random  search,  genetic 
algorithms,  Latin  hypercube  search,  or  orthogonal  arrays.  One  popular  strategy  is  to 
construct  and  optimize  a  surrogate  function.  A  surrogate  acts  as  an  approximation 
to  the  objective  function,  but  at  a  fraction  of  the  cost.  The  use  of  surrogates  often 
yields  significant  improvement  in  the  objective  function  value  early  in  the  iteration 
process.  For  further  details  on  surrogate  functions,  the  reader  is  referred  to  [35]. 

If  the  SEARCH  step  fails  to  provide  an  improved  mesh  point,  the  POLL  step  is 
invoked.  The  POLL  step  is  more  rigidly  defined  and  evaluates  the  neighboring  mesh 
points  for  the  current  iterate.  Use  of  positive  spanning  directions  in  the  construction 
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of  these  neighboring  points  provides  the  theoretical  basis  for  the  convergence  of  GPS. 
The  poll  set  at  iteration  k  can  be  expressed  as 


{xk  +  A kd  :  d  G  Bk},  (4.3) 

where  C  D  is  also  a  matrix  whose  columns  positively  span  Mn.  The  poll  set  is 
therefore  composed  of  mesh  points  neighboring  the  current  iterate  xk  in  the  directions 
of  the  columns  of  D&,  a  multiple  Ak  away  from  the  current  iterate. 

If  the  SEARCH  and  POLL  step  both  fail,  the  incumbent  solution  is  said  to  be  a 
mesh  local  optimizer  and  the  mesh  is  then  refined  by  setting  the  mesh  size  parameter 

Afc+1  =  0w»Ak,  (4.4) 

where  9  >  1  is  rational  and  wk  G  {w~  ,w~  +  1, . . . ,  —1}  for  some  w~ .  An  incumbent 
point  xk  is  replaced  by  xk+i  only  if  f(xk+i)  <  f(xk),  and  xk+\  is  termed  an  improved 
mesh  point.  If  an  improved  mesh  point  is  found  in  either  step,  then  the  mesh  is 
either  retained  or  coarsened  by  increasing  the  mesh  size  parameter  according  to 
equation  (4.4)  for  some  wk  G  {0, 1, . . . ,  iu+}.  It  follows  that  for  any  k  >  0  there 
exists  an  integer  rk  G  Z  such  that  Ak  =  6Tk  A0.  A  simplified  GPS  algorithm  for 
maximization  is  presented  in  Algorithm  1,  in  which  wk  =  0.  Algorithm  1  is  presented 
as  a  maximization  problem,  which  in  general  can  simply  be  considered  by  minimizing 
the  negative  value  of  an  objective  function  one  is  seeking  to  maximize. 

The  convergence  analysis  of  pattern  search  is  well-established  in  [21]  and  [130] 
and  requires  several  assumptions.  In  order  to  discuss  the  convergence  properties 
of  GPS,  several  assumptions  are  necessary.  First,  all  iterates  produced  by  GPS 
must  lie  in  a  compact  set  [47].  This  very  common  assumption  holds  as  long  as 
(xGfi:  f(x)  <  f(x0)}  is  compact.  Second,  if  the  matrix  G  =  I  (as  is  usually  the 
case),  then  the  constraint  matrix  A  must  be  rational.  The  final  necessary  assumption 
is  that  f(x o)  <  oo  for  xq  G  M™.  Torczon  [130]  proved  that,  under  the  assumptions 
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Let:  D  = 

-  {di,  d-2,  ■  ■  ■ 

,d\D\}  C  Mn  be  a  positive  spanning  set 

Given:  xq  G  Mn,  A0 

>  o,  e  >  i 

for  k  =  0, 

1, . . .  do 

if  3  de 

D  such  that  f(xk  +  Akd)  >  f(xk)  then 

Xk+1  < 

—  xk  +  A  kd 

A-k+1 

< 

1 

else 

Xk+ 1 

—  xk 

Afc+i 

<-  (l/0)Ak 

end  if 

end  for 

Figure  4.1  Algorithm  1:  A  basic  GPS  maximization  algorithm  for  unconstrained 
problems. 

discussed  above,  the  mesh  size  parameter  satisfies  lim  inf^+oc  A*,  =  0.  This  leads 
to  the  main  convergence  result  proved  in  [22]  and  restated  as  Theorem  4.1. 

Theorem  4.1  If  x  is  any  limit  of  a  refining  subsequence,  and  if  d  is  any  direction 
in  D  for  which  f  at  a  poll  step  was  evaluated  infinitely  often  in  the  subsequence, 
and  if  f  is  Lipschitz  near  x,  then  the  generalized  directional  derivative  of  f  at  x  in 
the  direction  d  is  nonnegative. 

Theorem  4.1  is  a  directional  result  only,  and  is  unsufficient  for  convergence  to  a 
stationary  point.  However,  since  the  optimization  problem  considered  in  this  thesis 
has  only  one  decision  variable,  D  =  {—1,1}  is  the  entire  set  of  normalized  positive 
spanning  directions  in  M  and  Theorem  4.1  applies  to  both  of  them.  Audet  [20]  proved 
convergence  to  a  Clarke  [46]  first-order  stationary  point  in  the  one- dimensional  case 
for  unconstrained  problems.  A  similar  argument  based  on  [14]  can  be  made  to  show 
convergence  under  mild  conditions  in  the  one-dimensional  case  to  a  local  minimizer. 

GPS  was  extended  in  [78]  and  [79]  to  problems  with  bound  and  linear  con¬ 
straints,  respectively.  To  handle  these  constraints  while  maintaining  convergence 
properties,  infeasible  points  are  discarded  without  being  evaluated,  and  search  di¬ 
rections  are  chosen  so  as  to  conform  to  the  geometry  of  the  nearby  constraint  bound- 
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aries.  In  one  dimension,  conforming  directions  are  moot,  since  there  are  only  two 
normalized  directions.  Extending  GPS  to  general  nonlinear  constraints  requires  aug¬ 
mentation  of  the  algorithm,  such  as  in  [80],  [24],  and  [23];  however,  in  one  direction, 
this  becomes  moot  for  the  same  reason  as  that  for  linear  constraints.  Thus,  in  one 
dimension,  it  is  sufficient  to  use  the  basic  GPS  algorithm  while  ignoring  infeasible 
points. 

The  NOMADm  optimization  software  [11,  13],  written  in  Matlab®  was  used 
to  implement  the  pattern  search  procedure  used  in  this  thesis.  NOMADm  is  specih- 
cally  designed  to  numerically  solve  nonlinear  and  mixed  variable  optimization  prob¬ 
lems  via  an  implementation  of  the  class  of  mesh  adaptive  direct  search  (MADS) 
algorithms.  GPS  is  a  subclass  of  MADS  [25],  in  which  poll  directions  are  restricted 
to  a  uniformly  bounded  finite  set.  In  fact,  in  one  dimension  MADS  and  GPS  essen¬ 
tially  coincide. 

The  mathematical  model  presented  in  chapter  3  must  be  converted  into  a 
MATLAB  function  that  produces  A(r)  for  the  NOMADm  optimizer.  This  is  ac¬ 
complished  by  first  directly  implementing  the  equations  of  chapter  3  into  a  single 
Matlab  function.  However,  naively  computing  the  equations  in  their  necessary  or¬ 
der  proved  too  computationally  expensive  and  resulted  in  excessive  run  times.  This 
motivated  the  need  to  investigate  a  number  of  computational  improvements  which 
are  discussed  in  section  4.3. 

Jy.2  Numerical  Computation  of  A(r) 

Recall  from  chapter  3  that  computation  of  A(r)  requires  values  for  p ,  E*  [iq] 
and  Ej  [Rh],  First,  a  direct  analytical  implementation  is  presented  that  simply  com¬ 
putes  the  equations  in  the  necessary  order  to  first  obtain  P,  which  leads  to  the 
vector  of  stationary  probabilities  p.  Subsequently,  the  conditional  expected  time 
to  the  first  failure  Ej  [R\]  and  first  replacement  E,  [Rh]  are  computed.  Direct  com- 
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putation  turns  out  to  be  extremely  computationally  expensive,  especially  when  the 
dimension  of  Q  or  the  value  of  A  is  large.  Thus,  several  computational  techniques 
were  developed  and  incorporated  in  an  attempt  to  improve  the  algorithm.  This  is 
followed  by  a  more  detailed  analysis  of  two  functions  which  account  for  the  majority 
of  the  computational  complexity:  the  matrix  exponential  and  Laplace  transform. 

In  order  to  compare  the  computational  complexity  of  various  algorithms  analy¬ 
sis,  techniques  have  been  developed  to  understand  complexity  growth  for  very  large 
input  values.  The  most  common  notation  used  in  complexity  analysis  is  referred  to 
as  big-oh  and  is  defined  as  follows  [121].  Given  any  two  functions  /  and  g,  then  f(n) 
is  0(g(n))  if  there  exist  positive  constants  C  and  k  such  that: 

f(n)  <  Cg(n )  whenever  n  >  k,  and  {n,  k}  C  Z. 

While  big-oh  notation  provides  insight  into  the  relative  rates  of  growth  between  al¬ 
gorithms,  it  does  not  describe  the  actual  characteristics  for  best  and  average  cases 
for  individual  run  times.  First,  many  algorithms  are  simply  too  complex  to  analyze 
mathematically  since  they  are  often  comprised  of  many  sub-functions  implemented 
in  different  languages.  Second,  it  is  difficult  to  qualitatively  account  for  memory 
management  and  the  interaction  between  hardware  architectures  and  specific  codes. 
For  example,  complexity  analysis  cannot  factor  in  the  effect  on  paging  as  virtual 
memory  usage  grows.  Third,  asymptotic  complexity  analysis  does  not  give  insight 
into  the  average  run  time.  Fourth,  asymptotic  complexity  analysis  does  not  give 
any  indication  of  the  run  time  or  efficiency  of  a  particular  algorithm,  only  an  un¬ 
derstanding  of  how  complexity  growth  varies  with  increasing  input.  For  example, 
one  algorithm  may  be  extremely  fast  and  another  cripplingly  slow  at  a  particular 
computation  but  have  the  same  algorithmic  order.  Because  of  the  disadvantages 
discussed  above,  only  computational  benchmarks  comparing  the  actual  run  times  of 
distinct  algorithms,  when  combined  with  asymptotic  complexity  analysis,  provide 
good  insight  into  the  relative  performance  between  algorithms. 
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Directly  computing  the  mathematical  model  presented  in  section  3.1  to  solve 
the  nonlinear  program  is  extremely  computationally  expensive.  The  full  algorithm  is 
listed  in  Algorithm  2.  It  is  comprised  of  five  different  stages,  each  with  a  particular 
computational  result.  The  first  stage  populates  the  transition  probability  matrix  P 
according  to  equation  (3.7).  The  complexity  of  this  stage  is  0(l2ryp(l))  if  0(p(l)) 
is  the  state-dependent  computational  burden  for  the  particular  matrix  exponential 
algorithm  needed  to  calculate  n This  stage  is  the  most  computationally  complex 
in  the  algorithm.  Within  the  third  nested  loop  there  are  two  Laplace  inversion  calls 
and  one  matrix  exponential  call  on  lines  6  and  7  of  Algorithm  2,  which  account 
for  the  primary  computational  expense  of  the  entire  algorithm.  For  this  reason, 
algorithms  for  these  two  function  calls  are  examined  carefully  in  section  4.4. 

The  complexity  analysis  above  reveals  sensitivity  to  the  state-size  when  build¬ 
ing  the  transition  probability  matrix  P.  For  each  element  of  P,  equation  (3.7) 
must  be  computed.  The  variable  7  defined  by  equation  (3.8)  represents  the  number 
of  inspections  that  can  occur  if  the  stochastically  degrading  component  reaches  its 
maximum  possible  lifetime.  The  maximum  lifetime  A  is  xjr\.  Therefore  if  A  is  much 
larger  than  r,  7  is  also  large.  Since  equation  (3.7)  involves  a  summation  of  7  terms, 
the  overall  computational  expense  is  very  sensitive  to  this  quantity. 

Stage  2  of  the  algorithm,  described  on  lines  13-16,  computes  the  stationary 
distribution  p  of  P.  This  is  simply  the  solution  to  the  equations, 

71  j  =  TTiPij  =  1 

i£S  i£S 

presented  in  [72],  These  steps  are  easily  implemented  and  do  not  incur  much  com¬ 
putational  cost. 

Stages  3  and  4  compute  equations  (3.13)  and  (3.14)  which  produce  E*  [F\] 
and  Ej  [i?i],  respectively.  Equation  (3.13)  is  implemented  on  lines  17  to  19  and 
equation  (3.14)  is  implemented  on  lines  20-26.  The  primary  driver  for  computational 
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Input:  r 

7  <-  \x/{rlT)~\ 

for  i  —  1  to  l  do 
for  k  —  1  to  l  do 

for  j  —  2  to  7  do 
7r  «—  expm(Qjr) 

Xj  <-  ^C_1  jr)|  -  £_1  ju_i(i 

end  for 
7r  <—  expm(Qr) 

P*,fc  <-  Y;X  +  (£_1  {M_1Gi(M,r)|j 

end  for 
end  for 
P<-P-I 

Set  all  elements  in  the  last  column  of  P  to  one 
Initialize  p  and  set  the  last  element  to  one 
p  <—  p  P-1 
for  i  —  1  to  l  do 


C  1  \  u  1Gi(u,  (j  -  l)r)  i  J  7Tit 


Set  the  ith  element  of  E  [ip]  to  £  1  yu  1ei  yuHo  —  Q  —  A(F d(u)  —  I)) 

end  for 

for  i  —  1  to  l  do 

for  j  =  1  to  7  —  1  do 

Sj  <-  £_1  {w_1Gi(-u,  jr)  j 

end  for 

Set  the  ith  element  of  E  [Ri\  to  r( 7  —  ^  S) 

end  for 

A  <-  (pEIFjD/fpEfflj]) 


Figure  4.2  Algorithm  2:  Direct  algorithm  for  evaluating  A(t). 


complexity  in  these  two  stages  is  the  Laplace  transform  inversion  algorithm  which 
is  examined  in  section  4.4.  The  complexity  of  stage  1  is  much  greater  than  either  of 
these  stages,  even  though  the  asymptotic  complexity  of  the  numerical  computation  to 
evaluate  equation  (3.14)  is  only  of  order  0(l( 7  — 1)).  The  final  stage  of  Algorithm  2  is 
to  bring  together  E,  E,  [iq]  and  p  to  compute  the  limiting  average  availability  A 
for  a  particular  r.  Equation  (3.4)  is  computed  on  line  26  without  much  computational 
expense. 

4-3  Computational  Enhancements 

Algorithm  2  was  presented  as  the  objective  function  for  NOMADm  and  several 
numerical  experiments  were  performed,  each  resulting  in  excessively  long  run  times, 
even  for  small  state  sizes.  This  motivated  the  need  for  a  more  efficient  implemen¬ 
tation  of  the  objective  function.  This  section  documents  several  strategies  taken  to 
reduce  the  computational  complexity  of  Algorithm  2,  since  this  should  reduce  overall 
run  times  for  the  optimization  procedure. 

The  first  strategy  was  to  focus  on  implementing  a  series  of  general  compu¬ 
tational  improvement  techniques.  Matlab  is  an  interpretive  language  with  more 
computational  overhead  than  an  associated  compiled  language.  Consequently,  there 
are  several  techniques  that  take  advantage  of  the  way  Matlab  executes  code  to 
vastly  improve  efficiency.  These  techniques  include  the  use  of  functions,  vectoriza- 
tion,  and  use  of  the  most  appropriate  operator. 

One  early  improvement  made  was  to  place  all  code  into  functions.  MATLAB 
code  executes  more  quickly  when  implemented  in  a  function  rather  than  a  script. 
Every  time  a  script  is  used,  the  entire  script  is  loaded  into  memory  and  evaluated 
one  line  at  a  time  with  a  large  amount  of  overhead  processing.  Functions,  on  the 
other  hand,  are  compiled  into  pseudo-code  and  loaded  into  memory  once,  with  pri- 
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vate  variables  accessible  only  to  the  current  function.  Therefore,  the  more  often  a 
procedure  is  called,  the  greater  the  advantage  of  using  a  function  instead  of  a  script. 

As  Matlab  is  designed  to  operate  on  matrices,  it  is  important  to  use  vector 
operations  whenever  possible.  Vectorization  can  increase  computational  speed  by 
several  orders  of  magnitude.  For  example,  while  for  loops  are  intuitive,  they  take 
much  longer  to  execute,  since  each  line  must  be  sent  separately  to  the  processor.  If  a 
vector  operation  is  used  instead,  Matlab  processes  this  operation  as  a  single  com¬ 
mand.  Vectorization  was  most  helpful  in  the  Laplace  transform  inversion  algorithm 
where  several  for  loops  were  replaced  by  different  vector  commands. 

An  example  of  replacing  an  inefficient  operator  with  one  more  suited  to  the 
computation  at  hand  was  to  replace  the  inv  command  on  line  16  of  Algorithm  2 
with  the  matrix  right  division  operator.  This  exchanges  the  full  LAPACK  inversion 
algorithm  for  a  custom  Matlab  procedure  that  exploits  the  structure  of  the  matrix 
to  find  the  most  efficient  method  of  solving  the  linear  system  (such  as  Gaussian 
elimination  with  partial  pivoting  if  the  matrix  is  square). 

The  second  strategy  used  to  increase  computational  speed  was  to  carefully 
analyze  locations  of  heavy  computational  cost  and  more  efficiently  implement  them. 
This  included  relocating  computations  from  the  more  frequently  called  sections  of 
Algorithm  2  to  a  pre-processing  block  that  is  only  executed  once  per  function  call. 
NOMADm  provides  for  the  use  of  a  parameter  hie  that  loads  at  the  start  of  a  run,  and 
each  value  that  requires  computation  only  once  was  moved  to  this  hie.  Since  memory 
retrieval  is  much  faster  than  computations,  even  seemingly  small  computations  were 
moved  into  this  hie. 

The  Matlab  profiler  measures  the  computational  cost  of  each  line  of  code. 
When  applied  to  the  MATLAB  implementation  of  Algorithm  2,  line  7  was  discov¬ 
ered  to  be  a  key  computational  cost  with  two  different  inversion  calls.  Since  the 
difference  between  the  two  calls  was  one  iteration,  this  line  was  changed  to  evaluate 
one  inversion  call  and  to  take  the  differences  all  in  one  step  outside  the  third  nested 
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loop  of  stage  one  in  Algorithm  2.  This  led  to  a  subsequent  analysis  showing  that, 
by  adjusting  the  indices  of  the  inner  loop,  one  expm  call  could  be  removed  from  the 
second  nested  loop  of  this  same  stage.  A  similar  analysis  was  conducted  throughout 
all  sub-functions. 

While  these  strategies  focused  on  making  Algorithm  2  run  more  quickly,  the 
third  strategy  focused  on  changing  the  algorithm  itself  to  improve  efficiency.  From 
MATLAB  profiler  measurements,  it  was  clear  that  most  computational  cost  was  com¬ 
ing,  as  expected,  from  the  third  nested  loop  of  stage  1.  Upon  observation,  the 
difference  between  each  function  call  and  the  next  quickly  approached  zero  and  de¬ 
creased  exponentially.  This  formed  the  central  idea  for  the  third  strategy,  to  truncate 
this  series  at  a  given  tolerance  level  e,  effectively  turning  Algorithm  2  from  a  direct 
method  to  an  iterative  one.  The  resulting  truncation  method  provided  savings  pro¬ 
portional  to  increasing  l  and  A.  The  truncation  method  can  dramatically  reduce  the 
number  of  terms  which  must  be  summed  by  approximating  the  (i,  fc)th  element  of 
P  according  to 

7 

Pi*  =  y^^i,k{nT)^j[Ri = nr\ 

n=  1 
N 

~  7TiAnr)  (' Gi(x >  nr )  -  Gi(x,  (n  -  l)r)) 

71=  1 

where  N  <  7  is  chosen  so  that  Gi(x,  Nr)  —  Gi(x,  (N  —  1)t)  <  e.  A  sensitivity  analysis 
was  conducted  and  the  value  e  =  10~4  was  chosen  in  order  to  improve  computational 
speed  to  the  maximum  degree  possible,  while  not  introducing  noticeable  error. 

4-4  Matrix  Exponential  and  Laplace  Transforms 

The  final  strategy  for  improving  computational  performance  was  to  analyze  the 
two  most  costly  functions:  the  matrix  exponential  and  the  inverse  Laplace  transform. 
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1:  Input:  r 

2:  Precalculate:  E  [Pi]  for  all  states 
3:  Set:  7  =  \x/{t1t)'\,  p 
4:  for  w  —  1  to  l  do 

5:  Clear  at  and  set  the  nth  element  of  ot  to  1 

6:  for  k  —  1  to  l  do 

7:  for  n  =  1  to  7  do 

8:  77  < —  Matrix  Exponential  of  Q t 

9:  Xn  Cr1  |n_1(5j(n,  nr) | 

P n  *  ^w,k 

11:  if  Xn  -  Xn_i  < e  then 

12:  Break 

13:  end  if 

14:  end  for 

15:  pw,k  y^2-o  [(Xr  ~  Xr_i)Vr\  +  V\Xi 

16:  end  for 

17:  end  for 

18:  P  ^  P-  I 

19:  Set  each  element  of  the  last  column  of  P  to  one 

20:  p  <—  p/P 

21:  for  i  —  1  to  l  do 

22:  for  n  =  1  to  (7  —  1)  do 

23:  Sn  =  £_1  |n_1Gi(n,nr)| 

24:  end  for 

25:  Set  the  ith  element  of  E  [Ri]  to  r( 7  —  ^  S) 

26:  end  for 

27:  Availability  * - (pE  [Pi])/(pE  [Pi]) 

28:  Cost  (p  (cd(E  [Ri\  -  E  [Pi])  +  cj|E  [Pi]  /rj  +  cRe))/(p E  [Pi])  -  B 
Figure  4.3  Algorithm  3:  Improved  algorithm  to  compute  A(r). 
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The  following  subsection  reviews  the  associated  numerical  methods  for  each  of  these 
functions. 

4.4.I  Computation  of  the  Matrix  Exponential 

In  order  to  numerically  calculate  the  availability  measures  discussed  in  section 
3.1,  it  is  necessary  to  obtain  the  transition  probability  matrix  for  {£n  :  n  >  0}, 
which  is  computed  according  to  equation  (3.7).  The  matrix  exponential  is  also 
a  necessary  operation  to  compute  a  component  lifetime,  per  equation  (3.3).  The 
matrix  exponential  is  defined  by 

eQ‘  =  E^-'  <«> 

k= 0 

In  order  to  compute  the  values  of  P  per  equation  (3.7),  equation  (4.5)  must 
be  numerically  computed  for  each  unique  combination  of  states  and  for  each  unique 
value  of  t  from  1  up  to  7.  The  limiting  complexity  of  matrix  exponential  evaluations 
is  0{l2 7)  or  potentially  varying  cubically  with  the  input.  Therefore,  for  environments 
with  even  moderate  state  sizes  (e.g.,  more  than  5)  and  moderately-sized  values  of 
7,  the  matrix  exponential  is  a  key  driver  of  computational  requirements.  For  this 
reason,  a  thorough  analysis  of  the  matrix  exponential  was  conducted. 

Numerous  techniques  for  computing  the  matrix  exponential  have  been  pre¬ 
sented  in  the  literature  [99,  133,  72,  128,  42,  92],  Moler  and  Van  Loan  [92]  in 
particular  surveyed  nineteen  different  methods.  A  common  method  shared  by  many 
solutions  uses  matrix-scaling  and  powering  to  first  form  Ordinary  differential 
equation  (ODE)  solvers,  can  compute  the  matrix  exponential  directly  without  ex¬ 
plicitly  forming  According  to  [135],  the  stability  of  the  matrix  exponential  is 
sensitive  to  perturbations  in  Q,  with  the  norm  of  Q  being  the  most  critical  measure 
of  stability.  According  to  [135],  if  Q  is  a  defective  matrix  (i.e.,  Q  does  not  have  a  full 
set  of  eigenvectors),  then  the  matrix  exponential  is  relatively  poorly  conditioned. 
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Several  methods  were  examined  for  computing  the  matrix  exponential,  each 
with  a  distinct  set  of  properties.  Specifically,  these  methods  include  a  simple  Tay¬ 
lor  series  expansion,  the  method  of  eigenvalues  and  eigenvectors,  the  uniformization 
method  and  the  Pade  approximation  via  repeated  scaling  and  squaring.  Computing 
the  matrix  exponential  using  these  four  different  methods  serves  two  primary  pur¬ 
poses.  First,  since  the  matrix  exponential  is  critical  to  the  overall  cost  of  computing 
A(t),  it  is  important  to  use  the  most  stable  and  efficient  method  available.  Second, 
by  using  four  different  methods,  it  is  possible  to  compare  results  and  mitigate  the 
risk  of  a  particular  method  causing  unacceptable  error  propagation  or  instability. 

Computation  of  the  matrix  exponential  via  Taylor  series  is  the  most  simple,  but 
often  slowest  and  least  accurate  method  studied.  According  to  the  classic  text  [57], 
it  is  only  useful  for  theoretical  comparisons.  Moler  and  Van  Loan  [92]  demonstrate 
that  this  implementation  is  subject  to  catastrophic  cancellation  as  truncation  errors 
are  greatly  magnified  as  the  factorial  of  the  denominator  increases.  Despite  these 
shortcomings,  the  simplicity  of  the  Taylor  series  method  facilitates  error  analysis 
and  provides  a  lower-bound  from  which  to  understand  other  methods.  A  simple 
expression  from  which  to  calculate  bounds  on  the  Taylor  series  truncation  error  is 
given  in  [81].  Furthermore,  the  Taylor  series  method  illustrates  the  classic  definition 
for  the  matrix  exponential  and  exposes  pitfalls  of  computational  round-off  error. 

The  theory  of  the  Taylor  series  method  is  straightforward.  From  elementary 
calculus,  for  any  matrix  T  and  scalar  t  one  may  compute  eTt  by  the  sum  of  powers, 


^  k\ 


(4.6) 


The  most  simple  way  to  deal  with  equation  (4.6)  is  to  add  terms  until  there  is  no 
machine-detectable  contribution  by  subsequent  terms.  The  marginal  contribution  of 
each  term  is  defined  by  the  matrix  F  and  the  cumulative  sum  of  the  expansion  is 
defined  by  the  matrix  E.  By  establishing  a  non-strict  boundary  and  using  the  matrix 
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F 

<-  I 

k 

<-  1 

while 

||F|| x  >  0  do 

E  <- 

-  E  +  F 

F  <- 

QF/k 

k  <- 

k  +  1 

end  while 

Figure  4.4  Algorithm  4:  Matrix  exponential  via  Taylor  series. 

1-norm  (i.e. ,  the  largest  column  sum),  the  condition  || F || x  >  0  stops  computation 
when  all  components  of  F  are  zero  to  machine  precision.  The  Math  Works  corporation 
provides  a  simple  MATLAB  implementation  of  the  Taylor  series  method  which  is 
summarized  in  Algorithm  4.  Results  of  the  Taylor  series  method  for  five  different 
example  problems  are  presented  in  chapter  5. 

Matrix  exponentiation  via  eigenvalues  and  eigenvectors  is  discussed  in  [72] 
and  [128].  It  is  most  stable  with  a  small  state  size  and  when  T t  has  distinct  eigen¬ 
values  [128].  Unlike  the  Taylor  series  and  uniformization  methods,  the  method  of 
eigenvalues  and  eigenvectors  can  handle  large  values  of  t,  as  is  often  the  case  for  re¬ 
liability  applications.  However,  a  small  matrix  size  is  preferred,  since  the  size  of  the 
matrix  significantly  impacts  the  cost  of  determining  eigenvectors  and  eigenvalues. 

The  method  of  eigenvalues  and  eigenvectors  is  built  on  the  relationship  that 
the  matrix  exponential  of  a  diagonal  matrix  is  the  diagonal  matrix  of  element  expo¬ 
nentials.  Consequently,  the  matrix  exponential  is  easily  computed  if  the  eigenvalues 
of  Q  are  known.  Assume  Q  €  Mnxn  has  n  distinct  eigenvalues  Uj,  j  =  1,2, .. . ,  n. 
Then  by  definition  [89],  if  Q  is  not  defective  (i.e.,  it  has  a  full  set  of  eigenvectors), 
then 

QjSj  "rs>  J  =  !,2, 
or 

Q  =  STS-1,  (4.7) 
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where  T  =  cliag  {u\ ,  z/2,  ■  ■  • ,  vn}  and  S  is  the  full  matrix  of  eigenvectors  in  where 
the  jth  column  is  the  right  eigenvector  Sj  corresponding  to  the  eigenvalue  Uj,  j  = 
1,2 ,  ...,n.  Exploiting  the  properties  of  diagonal  matrices,  equation  (4.7)  can  be 
written  more  generally  [72]  as, 


q  =  srKs 


fcc-l 


Accordingly,  equation  (4.5)  can  be  written  as, 


r  rtr 


eQ^SS-i  +  VS^S-1 

y  I 


=  s 


r= 1 
oo 


r  rtr 


r= 1 


i-l 


=  Ser‘S~‘. 


(4.8) 


(4.9) 


Equation  (4.9)  is  simple  to  calculate,  since  P  —  diag  {f[,  nj, . . . ,  id)}  ,  and  since  the 
matrix  exponential  of  a  diagonal  matrix  is  the  diagonal  matrix  of  element  exponen¬ 
tials,  it  follows  that 

ert  =  diag  {ePlt,  eV2\  . . . ,  eUnt}  .  (4.10) 

In  MATLAB,  the  eigenvalues  and  eigenvectors  can  be  easily  computed  using  the  built- 
in  function  eig,  which  uses  a  series  of  LAPACK  routines  to  compute  the  eigenvalues 
and  eigenvectors,  depending  on  the  properties  of  Q.  For  details,  see  the  LAPACK 
user’s  guide  [16]. 

The  main  advantage  of  the  method  of  eigenvalues  and  eigenvectors  [134]  is  that 
once  the  eigenvalues  and  eigenvectors  are  determined,  the  probability  distribution 
can  be  quickly  calculated  for  any  value  of  t  in  one  matrix-vector  multiplication. 
However,  in  this  research,  the  transient  solution  is  required  only  at  the  single  terminal 
point.  The  primary  weakness  for  this  method  is  that  Q  must  not  be  defective 
and  according  to  [128],  a  small  condition  number  cond(S)  of  S  (where  cond(S)  = 


4-19 


[V,  D]  <—  eig(Q) 

E  <—  V  diag  {exp(diag  {D})}  /V 

Figure  4.5  Algorithm  5:  Matrix  exponential  via  eigenvalues  and  eigenvectors. 

||S||  ||S-1||)  is  required  to  prevent  small  rounding  errors  from  compounding.  Moler 
and  Van  Loan  [92]  demonstrate  the  effect  of  large  condition  numbers  of  S  with 
numerical  examples.  A  MATLAB  implementation  of  this  method  is  summarized  in 
Algorithm  5.  Results  from  the  implementation  of  this  algorithm  are  presented  in 
chapter  5. 

The  uniformization  method  (also  called  Jensen’s  method)  is  thoroughly  dis¬ 
cussed  in  the  literature  [34,  72,  99,  128,  73,  85].  One  of  the  key  strengths  of  uni¬ 
formization  is  its  simple  implementation.  It  often  outperforms  other  methods  [128] 
and  performs  especially  well  when  Q  is  large  [73].  Furthermore,  uniformization  re¬ 
quires  only  modest  memory  allocation  and  provides  an  extremely  stable  and  efficient 
method  of  computing  the  matrix  exponential  [72],  The  primary  disadvantage  is  that 
large  values  of  ut  cause  considerable  roundoff  errors.  According  to  Kulkarni  [72],  if 
ut  >  250,  e~vt  is  potentially  less  than  the  smallest  computable  floating  point  number 
which  is  often  on  the  order  of  10~300. 

The  key  requirement  for  uniformization  is  that  the  diagonal  elements  of  Q  be 
bounded;  i.e. ,  if  rjjj  is  the  jth  diagonal  element  of  Q,  j  =  1,  2, . . . ,  n,  then 


\\Vjj\\  <  v  <  oo  (4-11) 

must  hold  for  each  state  j  and  for  some  v  >  0.  If  equation  (4.11)  holds,  Q  is  said  to 
be  uniformizable,  which  is  always  the  case  if  the  state  space  is  finite. 

Kulkarni  presents  the  following  procedure  in  [72],  First,  set  v  to 

v  =  max  {—rja}  ,  (4.12) 

\<i<i  L  1  v  ' 
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and  define  the  discretized  stochastic  transition  probability  matrix  Y  as 

T  =  I  +  iQ,  (4.13) 

which  produces  Q  =  z/(Y  —  I).  Since,  by  the  properties  of  exponentials, 

eQt  _  etvT-tvi  _  e~tv  e(tv) y 

the  unifomnization  equation  is  expressed  as 

OO  / 

T(f)=exp{Qt}  =  (4.14) 

k= 0 

where  the  inhnite  sum  in  equation  (4.14)  contains  the  Poisson  probability  mass 
function, 

Pr  {N(t)  =  k}  =  exp  (-id) . 
k\ 

A  distinct  advantage  of  uniformization  is  that  it  lends  great  control  over  trun¬ 
cation  error.  Stewart  [128]  presents  theoretical  bounds  for  the  truncation  error,  yet 
due  to  machine  round-off  errors,  the  computed  result  may  be  much  greater  than 
these  bounds.  However,  Stewart  [128]  states  that  roundoff  error  is  generally  not  a 
problem  since  numerical  operations  involve  no  negative  numbers  without  subtraction 
and  ||  Y||  <  1.  In  order  to  determine  at  which  value  M  to  truncate  the  sum,  so  that 
error  does  not  exceed  some  value  e,  the  Poisson  distribution  can  be  used  [72]  in  the 
equation, 

(4.15) 

k= 0 

The  uniformization  method,  shown  in  Algorithms  6  and  7,  was  implemented 
in  Matlab  as  presented  in  [99]  and  [128].  Algorithm  6  determines  the  value  of  M  in 
lines  4  through  9,  the  value  of  Y  according  to  equation  (4.13)  in  line  10,  and  the  value 
of  v  according  to  (4.12)  in  line  3.  Algorithm  6  is  called  only  once  per  calculation 
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of  A(t).  The  uniformization  equation  (4.14)  is  implemented  in  Algorithm  7  and  is 
called  with  every  matrix  exponential  call.  Further  code  was  implemented  on  the 
recommendation  of  [99]  to  condition  the  output  to  ensure  that  T  retains  properties 
of  a  stochastic  matrix.  Stewart  [128]  shows  that  the  complexity  of  the  uniformization 
method  is  0(n2)  and  requires  M(l  +  nz)  multiplications,  where  nz  is  the  number  of 
nonzero  elements  in  T. 


Input:  r,  Q,  £ 

Set:  K  =  0,  v  —  1  and  <7=1 
v  max(—  diag  {Q}) 

V  (i  -£)/(exp(-^r)) 
while  cr  <  r]  do 

K  <-  K  +  1 

v  (vut)  / K 

<7  < —  <7  V 

end  while 

Y  +  Q) _ 

Figure  4.6  Algorithm  6:  Determine  uniformization  parameters. 


Input: 

r,  v 

Y 

9  ) 

Set:  y  =  I 

for  k  = 

=  1  to  71  do 

y  <- 

(yYvr)/k 

7T 

n  +  y 

end  for 

E  exp  (—vt)  7t 


Figure  4.7  Algorithm  7:  Matrix  exponential  via  Uniformization. 


The  final  matrix  exponentiation  method  is  the  Pade  approximation  via  re¬ 
peated  scaling  and  squaring.  This  method  is  the  default  executed  by  the  Matlab 
expm  command  and  is  the  one  recommended  by  Moler  and  and  Van  Loan  [92],  The 
Pade  approximation,  without  a  method  to  control  roundoff  error,  suffers  the  same 
fate  as  the  Taylor  series  method  in  that  roundoff  error  worsens  as  ||Q||  t  increases  [27]. 
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This  error  can  be  mitigated  by  exploiting  a  property  of  the  exponential  function, 


exp(Q)  =  exp(Q/n)n. 


(4.16) 


A  value  of  n  can  be  chosen  to  be  a  power  of  two  for  which  exp(Q/n)  can  be  accu¬ 
rately  computed  [92],  From  this,  the  matrix  exp(Q/n)n  can  be  formed  by  repeated 
squaring.  Thus,  this  mechanism  scales  Q  by  a  power  of  two  to  reduce  the  norm 
to  order  one,  computes  a  Pade  approximation  to  the  matrix  exponential,  and  then 
repeatedly  squares  to  undo  the  effect  of  scaling. 

A  development  of  the  Pade  approximation  follows.  The  (a,  b )  Pade  approxi¬ 
mation  to  the  matrix  exponential  is,  by  definition  [17],  the  unique  rational  function, 

KMt)  s  (417) 
corresponding  to  constants  a  and  b.  The  (a,  b)  approximants  are  known  as  [61], 


Na,b(Qt)  =  Y, 

3=0 


(a  +  b  —  j)\a\ 
(a  +  b)\j\(a  —  j)\ 


(4.18) 


and 


b 


Da,b(Qt)  —  Y 

3=0 


(a  +  b  —  j)\b\ 
(a  +  b)\j\(b  —  j)\ 


(-Q  ty. 


(4.19) 


According  to  [128],  when  a  =  b  the  Pade  approximation  becomes  the  diagonal 
Pade  approximation  method,  in  which  case  equation  (4.17)  and  equation  (4.18) 
become 


R 


a, a 


Ng,a(  Qt) 
Naja(-Qt) 


(4.20) 


with, 


Na,a(Qt)  =  Y 

3=0 


(2  a-  j)\a\ 

(2  a)\j\{a-j)\ 


(Q  ty. 


(4.21) 
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The  diagonal  Pade  approximation  is  more  stable  than  other  choices,  since  in  Markov 
chain  problems,  the  eigenvalues  of  the  initial-value  problem  are  found  in  the  left 
half-plane,  causing  error  in  the  approximates  of  Ra,b{Qt)  [92],  If  a  >  b ,  cancelation 
problems  can  lead  to  large  roundoff  errors,  and  if  a  <  6,  Daib(Qt)  may  be  badly  con¬ 
ditioned  [92],  Furthermore,  a  higher-order  method  is  obtained  with  the  same  amount 
of  computation,  since  the  same  amount  of  work  is  needed  to  compute  Raa{ Qt),  and 
this  approximation  has  order  2 b  >  a  +  b. 

Equation  (4.21)  can  be  implemented  efficiently  as, 


1  +  2 


a/ 2-1 

(Qt)  ^  C2k+l(Qt) 

k= 0 


2k 


a/2  a/2—1 

^c2fc(Qt)2fc  -  (Qt)  ^2  C2fc+i(Qt) 


if  a  is  even 


2k 


Ran  =  < 


k= 0 


k= 0 


(a- 1)/2 

m  E  C2k  m 


(4.22) 


2k 


-1-2 


k= 0 


(a- 1)/2  (o— 1)/2 

^2  C2k+l(Qt)2k  -  ^2  C2k(Qi) 


if  a  is  odd 


2k 


k= 0 


k= 0 


where  c*  is  defined  by, 


c0  =  1 


Ci  —  C%— 1 


a  +  1  —  i 
i  (2 a  +  1  —  i)  ’ 


*  =  1,2,.... 


As  described  above,  a  major  disadvantage  of  using  the  Pade  method  is  that 
the  accuracy  of  the  approximation  decreases  with  distance  from  the  origin  and  thus 
degrades  with  the  increasing  size  of  ||Qt||2.  However,  if  for  a  and  b  are  sufficiently 
large,  or  if  the  eigenvalues  of  Qt  are  negative,  then  the  nonsingularity  of  Da  b(Qt )  is 
assured  [92],  Therefore,  Pade  approximation  can  be  used  if  || Qt ||  is  not  too  large. 

Moler  and  Van  Loan  [92]  perform  a  comparative  analysis  of  the  efficiency  of 
the  Pade  and  Taylor  approximants  and  analyze  their  relative  compatibility  with  the 
squaring  and  scaling  method.  They  find  that  the  combination  of  the  Pade  approx- 
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Require:  Scale  Q t  by  power  of  2  so  that  ||Ql||  <  1/2 

Input:  Q  t 

Set:  b  =  6,  a  =  True,  c  =  1/2 
{Compute  Pade  approximation  for  exp(Qt)} 

X^Q  t 
E  <—  I  +  cQt 
D  <—  I  —  cQt 
for  k  =  2  to  b  do 

c  <—  c(b  —  k  +  l)/k(2b  —  k  +  1) 

X  <-  Q  tX 
ip  <-  cX 
E  <—  E  +  ^ 
if  a  is  True  then 
D  <—  D  +  ip 
else 

D  ^  D  -ip 

end  if 

Set  a  to  the  Boolean  conjugate  of  its  current  value 

end  for 
E  <-  D/E 

{Undo  scaling  by  repeated  squaring} 
for  k  =  1  to  s  do 

E  <-  E2 

end  for 

Figure  4.8  Algorithm  8:  Pade  Approximation  via  squaring  and  scaling. 

imation  with  repeated  squaring  and  scaling  is  the  most  stable  and  computationally 
efficient  method.  Moreover,  the  Matlab  expm  command  employs  the  Pade  approx¬ 
imation  via  repeated  squaring  and  scaling.  Matlab’s  implementation  cites  [92]  and 
[59],  with  guidance  from  Ward  [140],  as  is  detailed  in  Appendix  A.  Their  procedure 
is  demonstrated  in  Algorithm  8. 


4-4-2  Numerical  inversion  of  Laplace  transforms 

In  order  to  compute  the  numerical  inversion  of  the  Laplace  transforms  in  equa¬ 
tion  (3.15),  a  method  provided  by  Abate  and  Whitt  in  [5]  was  chosen.  In  [5]  two 
inversion  algorithms  are  presented,  both  variants  of  the  Fourier-series  method.  The 
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Fourier-series  method  is  the  term  used  to  describe  the  process  of  implementing  the 
Fourier  series  of  an  associated  periodic  function  to  numerically  integrate  via  trape¬ 
zoidal  integration.  The  error  for  this  algorithm  is  bounded  only  partially  by  using 
the  Poisson  summation  formula  to  identify  the  discretization  error  associated  with 
the  alternating  series  in  the  algorithm  provided  by  the  trapezoidal  rule. 

The  inverse  Laplace  transform  is  defined  as, 


£_1 


pa+ioo 


es  f(s)  ds , 


(4.23) 


where  f(s)  is  the  Laplace  transform  of  the  function  f(t),  and  i  =  \f—l.  If  f(t)  is 
a  real-valued  function,  the  right-hand  side  of  equation  (4.23)  can  be  expressed  per 
[52]  as 
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(4.24) 


where  Re[-]  denotes  the  real  part  of  a  complex  number. 


In  practice,  explicitly  computing  equation  (4.24)  is  difficult  and  often  not  pos¬ 
sible,  requiring  the  use  of  numerical  methods.  There  is  a  significant  body  of  work  in 
the  literature  on  various  techniques  to  numerically  evaluate  equation  (4.24)  such  as 
[132,  8,  7,  6,  9,  5,  58,  3,  4,  2,  141,  1],  Most  of  these  articles  apply  the  finite  fourier 
cosine  transform,  which  relates  directly  to  the  Laplace  transform.  Dubner  and  Abate 
introduced  an  approximation  formula  in  their  seminal  paper  [52],  using  this  method 
to  approximate  equation  (4.24)  to  any  desired  accuracy  according  to  the  following 
formula, 
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where  gn(t)  is  one  of  an  infinite  number  of  even  periodic  functions,  each  with  pe¬ 
riod  2 T.  Equation  (4.25)  is  actually  the  the  trapezoidal  rule  [38]  applied  to  equa¬ 
tion  (4.24),  which  is  the  foundation  upon  which  most  numerical  Laplace  transform 
inversions  are  implemented. 

However,  numerous  other  methods  have  been  implemented  both  by  extending 
the  method  of  [52]  into  multiple  dimensions  and  attempting  to  use  other  methods 
such  as  Laguerre  functions  [141].  However,  methods  employing  Laguerre  functions 
perform  poorly  on  problems  with  large  values  of  t,  and  further  problems  are  intro¬ 
duced  by  their  non-geometric  convergence  [120]. 

In  this  thesis,  Laplace  transforms  are  inverted  using  an  algorithm  of  Abate  and 
Whitt  [5],  which  uses  a  variant  of  the  Fourier-series  method  specifically  tailored  to 
Laplace  transforms  of  probability  distributions.  Abate  and  Whitt  [5]  present  two 
distinct  methods  to  invert  equation  (4.24).  The  first  uses  the  Bromwich  integral 
[116],  the  Poisson  summation  formula,  and  the  Euler  summation,  while  the  second 
uses  the  Post-Widder  formula,  the  Poisson  summation  formula,  and  the  Stchfest 
enhancement  [2] . 

Equation  (4.24)  can  be  evaluated  by  the  trapezoidal  rule  with  step  size  h  = 
7t/2x  and  a  =  A/2x,  yielding  the  following  series, 


eA /2 

h(x)  =  ~^rRelf} 


(4.26) 


Since  equation  (4.26)  involves  an  infinite  sum,  Abate  and  Whitt  [5]  chose  the  Euler 
summation  acceleration  technique  due  to  its  simplicity  and  adequate  computational 
efficiency.  Defining  sn(x)  as  the  approximation  to  /),(x)  in  equation  (4.26)  via  Euler 
approximation  techniques,  sn(x)  becomes 
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where  dk(x)  is  defined  as 


=  Re[/]  (^^)  •  (4.28) 

Since  the  Euler  summation  is  applied  to  m  terms  after  an  initial  n,  the  full  Euler 
sum  approximation  is 

m 

E(m,  n,  x )  = 

fc= o 

and  Abate  and  Whitt  [5]  recommend  the  values  m  —  11  and  n  =  15.  Equations 
(4.27),  (4.28)  and  (4.29)  were  implemented  in  the  MATLAB  function  invt.m,  given 
in  Appendix  C. 

Because  the  Laplace  transform  inversion  algorithm  must  be  called  with  each 
function  call  in  the  main  availability  computation  (Algorithms  2  and  3),  it  is  critical 
to  minimize  the  computational  time  of  the  inversion  algorithm.  The  inputs  to  the 
inversion  procedure  are:  the  function  to  be  evaluated  (in  this  case,  the  component 
lifetime  distribution  function  Gi(x,t)),  the  initial  probability  distribution  vector  e, 
and  the  current  inter-inspection  interval  r.  Since  the  computational  complexity  of 
this  algorithm  does  not  vary  with  the  input  size,  it  is  0(1). 

Even  with  an  asymptotic  complexity  of  0(1),  (the  algorithm  was  implemented 
in  UBASIC  by  Abate  and  Whitt  in  [5])  it  contains  over  140  lines  of  code.  In  order 
to  reduce  this  computational  expense,  quantities  that  only  needed  to  be  computed 
once  were  pre-calculated  and  placed  outside  of  Algorithm  9.  In  [5],  Abate  and  Whitt 
recommended  values  for  the  constants  for  the  Euler  summation  acting  as  an  average 
of  the  last  m  partial  sums  of  a  binomial  probability  distribution  with  parameters  m 
and  p  =  0.5.  They  applied  the  Euler  summation  to  m  terms  after  an  initial  n  of 
equation  (4.29)  and  recommended  values  of  m  =  11  and  of  n  =  16.  For  notational 
simplicity  and  computational  savings,  they  set  the  quantity  211  in  equation  (4.29)  to 


m 


Sfi+k  (x)  i 


(4.29) 
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Input:  EG,  zO,  nr 

Load:  a,  u,  c,  fi,  h,  v 

G{x,a,nT,f{x))/ 2 

for  k  =  1  to  v  do 

<t  <r  +  (-l)fc  G(x  +  khi,  ol,  nr,  f(x  +  khi)) 

end  for 

Store  c;  as  the  first  element  of  a 

for  k  —  1  to  12  do 
n  v  +  k 

a(k  +  1)  <—  a(k)  +  (—1)”  G  (x  +  nhi,  ol,  nr,  f(x  +  nhi)) 

end  for 

Output  u(c  ■  o')/ ii 

Figure  4.9  Algorithm  9:  Laplace  transform  inversion  algorithm. 

Ii~l  and  A  =  81n(10).  The  quantity,  ( eAk2)/x  in  equation  (4.27)  is  assigned  to  the 
variable  u  and  the  quantity  x  —  A/{2x)  and  h  =  n/x. 

In  this  chapter,  typical  optimization  strategies  were  surveyed.  Generalized 
pattern  search  was  presented  as  a  suitable  optimization  algorithm,  and  a  specific 
Matlab  implementation  was  discussed.  A  direct  and  improved  algorithm  for  com¬ 
puting  A{r)  was  presented  with  a  discussion  of  three  different  computational  im¬ 
provement  strategies  applied  to  the  computation  of  A[r),  with  specific  emphasis  on 
the  matrix  exponential  and  numerical  Laplace  inversion. 

The  next  chapter  illustrates  the  implementation  of  the  solution  procedures 
presented  here.  After  a  description  of  the  overall  experiment,  the  pattern  search 
methodology  used  to  maximize  A(t)  will  be  discussed.  Also,  the  results  from  the 
matrix  exponential  study  will  be  presented.  Five  different  cases  will  be  considered, 
each  with  a  distinct  set  of  parameters.  Primarily  characterized  by  the  number  of 
environmental  states,  the  cases  were  designed  to  test  the  broad  applicability  of  each 
chosen  methodology.  A  particular  emphasis  is  placed  upon  the  computational  im¬ 
provements  obtained  in  each  case  by  the  use  of  the  truncation  method. 
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5.  Numerical  Results 


This  chapter  summarizes  and  discusses  the  results  of  a  series  of  numerical 
experiments  involving  the  maximization  of  A  by  selection  of  an  appropriate  inter¬ 
inspection  interval,  r.  All  numerical  computations  were  performed  in  MATLAB.  The 
NOMADm  implementation  of  the  GPS  algorithm  was  used  to  solve  for  the  approx¬ 
imate  optimum  inter-inspection  interval.  Results  are  presented  for  five  scenarios 
involving  degradation  from  environment-dependent  linear  wear  and  random  shocks. 

Two  experiments  were  designed  to  illustrate  the  optimization  problem  and  to 
evaluate  particular  computational  improvements.  The  first  experiment  evaluated 
four  different  implementations  of  the  matrix  exponential  algorithm  in  order  to  de¬ 
termine  the  most  efficient  and  stable  implementation.  The  meaning  of  these  terms 
in  the  context  of  this  research  will  be  discussed  below.  The  second  experiment 
examined  the  performance  gains  provided  by  the  truncation  method  discussed  in 
chapter  4.  The  truncation  method  used  only  the  most  suitable  matrix  exponential 
algorithm  resulting  from  the  previous  experiments. 

In  numerical  analysis,  the  stability  of  an  algorithm  refers  to  its  accuracy.  A 
given  algorithm  is  numerically  stable  if  it  produces  a  good  approximation  to  the  true 
solution  [77].  Since  solutions  produced  in  this  thesis  are  numerical  approximations, 
stability  will  refer  to  the  particular  capability  of  an  algorithm  to  arrive  at  a  solution 
consistently  produced  by  algorithms  of  known  stability.  Often  a  given  quantity  can 
be  numerically  computed  in  several  different  ways,  all  of  which  are  algebraically 
equivalent,  but  in  practice  yield  different  results  because  one  method  is  more  stable 
than  another.  An  algorithm  will  be  termed  unstable  if  it  does  not  return  an  answer, 
due  to  a  undesired  halt  in  code  execution.  If  one  algorithm  is  more  efficient  than 
another,  it  more  effectively  uses  computational  resources  (such  as  time  and  storage) 
in  order  to  perform  a  given  computation  [143] .  In  this  research,  efficiency  is  measured 
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by  a  combination  of  algorithmic  complexity  and  the  elapsed  run  time  required  to 
execute  an  algorithm. 

Both  experiments  were  conducted  in  the  MATLAB  computing  environment  on 
a  Silicon  Graphics  (SGI)  Origin  3900  maintained  by  the  Air  Force  Aeronautical 
Systems  Center.  The  Origin  3900  is  a  large-scale  server  running  the  custom  SGI 
unix-based  operating  system  IRIX  with  over  2048  processors  providing  a  peak  com¬ 
putational  power  of  2.9  x  1012  floating  point  operations  per  second  (FLOPS).  Each 
processor  is  a  SGI  MIPS  R16000  with  a  clock  speed  of  700  MHz  and  one  gigabyte 
of  dedicated  random  access  memory  (RAM).  A  user  is  given  exclusive  access  to  each 
allocated  processor  and  its  associated  RAM  for  computational  use.  This  is  useful 
when  comparing  the  computational  time  required  by  individual  runs,  as  there  are 
no  fluctuating  background  processes  interfering  with  benchmarks,  as  is  the  case  on 
a  personal  computer.  The  relative  machine  accuracy  for  floating  point  numbers  on 
the  Origin  3900  using  Matlab  is  2.2204E— 16.  The  primary  benchmark  to  measure 
a  given  algorithm’s  performance  was  the  clock-time  elapsed  during  each  individual 
run,  denoted  by  p. 

Five  different  cases  were  constructed,  each  with  a  distinct  combination  of  char¬ 
acteristic  parameters  to  fully  test  the  algorithms  presented  over  a  wide  range  of 
scenarios.  These  parameters  are:  the  generator  matrix  for  the  environment  with 
an  associated  linear  wear  rate  for  each  state,  the  shock  distribution  and  Poisson 
arrival  rate,  along  with  the  critical  damage  threshold  for  the  system.  Parameters 
were  chosen  in  order  to  demonstrate  the  optimization  algorithms  over  a  wide  range 
of  possibilities  and  are  presented  in  the  following  sections  in  order  of  the  associated 
dimension  of  the  random  environment. 

In  the  seven-,  ten-  and  twenty-state  cases,  the  transition  rates  present  in  the 
generator  matrix  Q  were  randomly  generated  using  the  MATLAB  rand  function. 
To  create  a  valid  infinitesimal  generator  matrix,  a  Matlab  script  set  the  diagonal 
elements  to  the  negative  of  the  sum  of  the  remaining  elements  in  the  row.  The 
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Matlab  rand  function  uses  a  modified  version  of  Marsaglia’s  [86]  subtract  with 
borrow  algorithm  and  can  theoretically  generate  over  21492  values  before  repeating 
itself.  Other  parameters  were  chosen  in  order  to  cover  a  range  of  distributions  and 
wear  rates  along  with  varying  values  of  the  maximum  lifetime  A  =  xjr\.  This  ratio  is 
of  particular  importance,  since  computational  requirements  are  extremely  sensitive 
to  A  when  combined  with  small  values  of  r. 

In  order  to  control  and  initiate  batch  runs,  a  series  of  scripts  in  the  Perl  com¬ 
puter  language  were  produced  to  prepare  the  hie  structure  for  each  run  and  interface 
with  the  MSRC  load  sharing  facility  (LSF)  to  submit  jobs.  Batch  runs  generating 
replications  were  necessary,  since  the  elapsed  clock-time  provides  the  computational 
basis  of  comparison  between  any  two  algorithms,  and  there  is  an  inherent  noise  due 
to  the  physical  characteristics  of  an  individual  processor  that  can  be  approximated 
as  random  noise  present  in  computational  run  times  [60].  The  effect  of  noise  becomes 
more  pronounced  as  the  number  of  different  processors  used  in  the  experiment  in¬ 
creases.  The  noise  present  in  computational  run  times  on  large-scale  servers  is  due 
to  interactions  from  an  number  of  processes  [60].  Most  prominent  among  these  are 
access  to  level-3  cache  and  differences  in  the  characteristics  of  the  various  processors 
performing  runs  in  the  server.  Since  a  large-scale  server  contains  a  number  of  dif¬ 
ferent  processors,  each  with  a  unique  actual  processing  speed,  there  is  an  inherent 
noise  in  results,  even  when  processors  are  dedicated  to  a  particular  user  and  all  with 
equal  processing  speeds. 

For  these  reasons,  a  statistical  analysis  was  conducted  to  determine  the  number 
of  runs  required,  and  the  main  batch  generation  script  generated  30  jobs  for  each 
combination  of  method  and  environmental  scenario.  This  number  was  chosen  in  or¬ 
der  to  provide  a  sufficiently  small  confidence  interval  for  the  computational  run  time 
while  still  working  within  the  budget  available  for  computing  resources.  After  com¬ 
pletion  of  all  experiments,  a  shell  script  collected  all  data,  which  were  then  processed 
and  analyzed  on  a  personal  computer  using  a  custom  Matlab  post-processing  script. 
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The  first  experiment  evaluated  four  different  methods  to  compute  the  matrix 
exponential  as  discussed  in  chapter  4.  These  methods  are  the  Taylor  series,  the 
method  of  eigenvalues  and  eigenvectors  (EE),  the  Uniformization  method  and  the 
Pade  approximation  via  repeated  squaring  and  scaling  (PSS).  Each  scenario  was 
examined  and  evaluated  for  stability  and  computational  efficiency.  In  order  to  de¬ 
termine  the  relative  efficiency  of  the  various  methods,  a  series  of  heteroscedastic, 
two-tailed  f-tests  were  conducted  to  test  the  null  hypothesis  that  there  is  no  statisti¬ 
cal  difference  between  the  run  times  of  any  two  given  methods.  Conhdence  intervals 
for  p  were  computed  at  a  0.05  level  of  significance.  After  completion  of  the  analysis, 
the  most  stable  and  computationally  efficient  matrix  exponentiation  method  was 
found  to  be  the  Pade  approximation  via  repeated  squaring  and  scaling. 

The  second  experiment  used  the  results  of  the  first  to  specifically  consider  the 
computational  advantage  provided  by  the  truncation  method.  As  in  the  first  exper¬ 
iment,  30  runs  were  conducted  for  the  baseline  and  truncation  method  for  all  five 
cases.  The  baseline  performed  the  full,  direct  computation,  while  the  truncation 
method  trimmed  all  successive  terms  smaller  than  0.0001.  Similarly,  a  heteroscedas¬ 
tic  two-tailed  f-test  was  conducted  to  ensure  any  improvement  in  computational  run 
time  due  to  truncation  was  statistically  significant. 

5.1  NOMADm  Configuration 

This  section  discusses  the  particular  configuration  of  NOMADm  used  to  pro¬ 
duce  the  results  of  the  numerical  experiments  presented  in  the  following  sections. 
The  overall  optimization  problem  was  set  up  to  work  with  a  main  function  hie, 
named  availabilityCalc,  which  evaluates  an  iterate  r  and  returns  the  correspond¬ 
ing  objective  A(t)  and  cost  function  values.  A  hie  named  availabilityCalc_xO 
simply  returns  the  initial  point  which  is  the  center  of  the  maximum  lifetime  (A)  of 
the  component.  The  feasible  region  (with  respect  to  bound  constraints)  for  the  deci¬ 
sion  variable  is  defined  in  a  hie  named  availabilityCalcJDmega.  In  order  to  reduce 
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computational  complexity,  all  calculations  requiring  only  one  evaluation  were  placed 
in  a  separate  parameter  file  availabilityCalc_Param,  which  returns  a  Matlab 
structure  whose  fields  contain  all  parameters  used  to  characterize  a  particular  sce¬ 
nario,  plus  any  other  calculations  that  may  be  computed  in  advance.  The  use  of  this 
file  had  a  dramatic  effect  on  decreasing  the  overall  run  time. 

The  one-dimensional  poll  directions  for  this  problem  were  chosen  as  {—1,1}, 
and  the  poll  set  was  evaluated  only  until  an  improvement  is  found,  as  opposed  to 
evaluating  the  entire  set.  The  initial  starting  point  To  was  chosen  to  be  the  midpoint 
of  the  feasible  region  (0,  A] ,  which  seems  reasonable  because  little  information  is 
available  concerning  the  structure  of  the  optimization  problem.  While  most  opera¬ 
tionally  important  values  of  r  are  small  relative  to  the  maximum  lifetime,  starting 
at  a  very  small  value  of  r  significantly  increases  computational  cost  unnecessarily. 
This  is  because  the  component  maximum  lifetime  A  remains  fixed,  and  small  values 
of  r  dramatically  increase  7  and,  consequently,  the  computational  time.  Since  the 
pattern  search  method  works  by  forming  a  mesh  which  is  reduced  in  size  as  the 
algorithm  progresses,  it  is  much  better  to  approach  the  optimal  solution  and  take 
large  steps  in  the  potentially  wrong  direction  where  r  values  are  larger. 

As  the  MADS  optimizer  implemented  by  NOMADm  is  an  iterative  method, 
termination  was  set  to  occur  when  the  mesh  size  shrunk  below  a  tolerance  of  10-4, 
based  upon  some  preliminary  empirical  studies,  or  when  the  number  of  function 
evaluations  exceeded  50,000. 

NOMADm’s  speed  of  convergence  is  sensitive  to  the  settings  governing  mesh 
control.  The  initial  mesh  size  was  set  to  A0  =  A/4  as  a  compromise  between  a  step 
sufficiently  large  to  quickly  approach  the  optimal  r,  but  not  so  large  as  to  overshoot 
and  end  with  a  small  mesh  in  the  computationally  expensive  region  where  r  is  small. 
The  mesh  refining  factor  was  set  to  0.5,  and  no  mesh  coarsening  was  employed. 

The  shock  distribution  parameters,  shock  arrival  rate,  and  critical  damage 
threshold  for  each  run  are  presented  in  exact  representations  in  Table  5.1.  Cost 
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Table  5.1  Key  run  parameters. 


States 

Shock  Distribution 

Shock  Mean 

A 

X 

min(r) 

A 

2 

Exp  (4) 

0.25 

0.5 

1 

0.25 

4 

5 

Er(0.2,8) 

40 

0.25 

100 

1 

100 

7 

U(0,10) 

5 

1 

50 

2.5 

20 

10 

Er(0.5,4) 

8 

1 

60 

0.4 

150 

20 

Gamma(2,4) 

8 

1 

100 

0.94 

106 

parameters  were  held  constant  in  both  experiments  and  in  all  five  runs.  The  cost  of 
downtime  per  unit  time  was  set  to  0.5  units  per  unit  time,  the  cost  per-inspection  was 
set  to  1  unit  and  the  cost  of  replacement  was  set  to  a  fixed  5  units.  The  generator 
matrices,  linear  wear  rates  and  Laplace-Stieltjes  transforms  for  the  shock  damage 
magnitude  are  presented  individually  in  the  sections  that  follow. 


5.2  2- State  Case 


The  two-state  case,  adapted  from  [67],  models  a  system  that  transitions  be¬ 
tween  two  environment  states  S  =  {1,2},  each  with  a  distinct  rate  of  linear  wear. 
Environmental  transitions  occur  according  to  the  infinitesimal  generator  matrix, 


Q 


-25/3  25/3 

25/3  -25/3 


and  the  following  linear  wear  rates, 


r  = 


13/12  1/4 


The  random  damage  magnitude  due  to  shocks  is  assumed  to  be  distributed  expo¬ 
nentially  with  rate  4,  and  the  Laplace-Stieltjes  transform  of  this  distribution  is  given 
by 


4.0 

4.0  +  u' 
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Figure  5.1  Availability  (-)  and  cost  (- -)  for  the  2-state  case. 


In  this  example,  shocks  arrive  according  to  a  Poisson  process  with  rate  para¬ 
meter  A  =  0.5.  The  critical  damage  threshold  level,  x,  was  set  to  1.0.  From  the  wear 
rates  presented  above,  the  minimum  wear  rate  is  0.25  and  the  damage  threshold  is 
1.0  producing  a  maximum  lifetime  of  4  time  units.  The  budget  is  set  to  35.  A  plot 
of  the  limiting  average  availability  and  cost  over  several  values  of  r  is  presented  in 
Figure  5.1. 

The  results  from  the  matrix  exponentiation  experiment  are  displayed  in  Ta¬ 
ble  5.2.  All  four  matrix  exponentiation  methods  produce  the  approximate  optimal 
inter-inspection  time  of  r*  =  1.69,  with  the  Taylor  series  method  slightly  out  of 
agreement  with  the  others.  The  Taylor  series  method  did  result  in  24  function  calls 
to  reach  an  approximate  optimal  solution  while  the  other  three  only  required  23 
which  suggests  instability  with  the  Taylor  series  method  as  discussed  in  chapter  4. 
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Table  5.3  displays  the  p- values  of  the  methods  relative  to  each  other.  From 
Table  5.2,  we  see  that  the  EE  method  requires  the  smallest  run  time,  followed  by  the 
Uniformization  method.  This  must  be  stated  cautiously  since  the  EE  method  uses 
a  quickly  executing  compiled  LAPACK  function  to  compute  the  eigenvectors  and 
eigenvalues.  According  to  Tabic  5.3  the  EE  method’s  computational  advantage  is 
statistically  significant  at  a  high  level  of  significance  (at  least  on  the  order  of  1CT10). 
Also  notable,  is  that  the  Taylor  series  is  the  statistically  the  slowest  method. 


Table  5.2  Two-state  matrix  exponentiation  results. 


Uniformization 

EE 

Taylor 

PSS 

T* 

1.6971 

1.6942 

1.6877 

1.6942 

A  (r*) 

0.68989 

0.68980 

0.68957 

0.68980 

P 

12.065 

11.960 

13.613 

12.094 

°p 

0.048449 

0.048628 

0.032390 

0.081434 

Table  5.3  Two-state  matrix  exponential  p- value  comparison. 


EE  Taylor  PSS 

Uniformization 

EE 

Taylor 

1.73E-11  4.96E-68  1.14E-10 
2.75E-69  7.13E-10 
9.93E-47 

For  the  2-state  case,  there  is  no  statistical  difference  between  the  computational 
run  times  when  contrasting  the  baseline  and  truncation  method  with  an  associate 
p- value  from  the  f-test  of  0.15.  Both  methods  found  the  approximate  optimal  r  in 
around  12  seconds.  This  is  expected,  since  A  and  the  dimension  of  Q  are  both  small, 
the  gain  provided  by  the  truncation  method  is  not  significant,  ft  is  also  important  to 
note  that  the  only  measure  showing  a  measurable  difference  among  the  thirty  runs 
is  the  computational  run  time.  There  is  no  detectible  variance  between  the  values 
of  r*  or  A(t*)  in  this  or  any  of  the  following  cases.  The  mean  run  time  (p)  for 
the  baseline  method  was  12.11  seconds  in  the  interval  from  12.09  to  12.12  seconds. 
For  the  truncation  method  p  was  12.09  seconds  in  the  interval  from  12.07  to  12.11 
seconds. 
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Table  5.4  Two-state  truncation  method  and  baseline  comparison. 


Full  Sum 

Truncated  Sum 

T* 

A(C) 

p 

1.694213867 

0.689801977 

12.10796217 

0.042446646 

1.694213867 

0.689801977 

12.09131087 

0.046830392 

5.3  5- State  Case 

In  the  5-state  case,  which  is  also  adapted  from  [67],  shocks  arrive  according  to 
a  Poisson  process  with  a  rate  parameter,  A  =  0.25.  The  critical  damage  threshold 
level  is  set  at  x  =  100,  and  the  overall  budget  for  inspections,  down-time,  and 
replacements  is  2  units.  The  generator  matrix  is, 


Q  = 


-0.500 

0.125 

0.125 

0.125 

0.125 

0.400 

-2.000 

0.400 

0.600 

0.600 

0.025 

0.025 

-0.100 

0.025 

0.025 

0.050 

0.050 

0.050 

-0.200 

0.050 

1.500 

1.000 

1.000 

1.500 

-5.000 

The  condition  number  of  Q  is  on  the  order  of  101'  which  could  cause  stability  prob¬ 
lems  with  the  matrix  exponential  calculations  using  either  the  Taylor  series  or  the 
EE  method.  The  vector  of  linear  wear  rates  corresponding  to  each  state  is 

r  =  [  1  2  3  4  10  ]• 


Since  the  minimum  wear  rate  is  one  and  x  =  100  the  system  lifetime  of  100  is  much 
longer  than  in  the  two-state  case.  This  causes  a  larger  computational  burden,  due 
to  the  much  larger  value  of  7  when  r  is  small.  The  shock  magnitudes  are  Erlang 
distributed  with  parameters  0.2  and  8,  such  that  the  corresponding  Laplaee-Stieltjes 
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Figure  5.2  Availability  (-)  and  cost  (--)  for  the  5-state  case. 


transform  of  the  shock  distribution  is 


Fy{u) 


(—) 

V0.2  +  uj 


A  plot  of  the  limiting  average  availability  and  cost  over  time  for  this  case  is  presented 
in  Figure  5.2. 


Table  5.5  Five-state  matrix  exponentiation  results. 


Uniformization 

EE 

Taylor 

PSS 

T* 

5.7961 

5.7961 

50.154 

5.7961 

A(r*) 

0.75477 

0.75580 

0.20703 

0.75580 

P 

867.64 

852.74 

85.81 

865.33 

oP 

4.34 

5.41 

0.672 

10.2 

Table  5.5  displays  the  five-state  case  matrix  exponential  results.  It  is  imme¬ 
diately  clear  that  the  Taylor  series  method,  while  much  faster,  is  producing  vastly 
different  results  which  implies  instability  and  consequently  the  Taylor  series  method 
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Table  5.6  5-state  p- values  for  matrix  exponentiation  performance  tests. 


EE  Taylor  PSS 

Uniformization 

EE 

Taylor 

1.1E-16  7.0E-70  0.26 

8.4E-66  3.4E-7 

7.7E-57 

cannot  be  considered  a  viable  option  for  this  case.  The  three  remaining  options 
required  73  function  evaluations  compared  to  the  27  required  by  the  Taylor  series 
method. 

The  p- values  resulting  from  comparative  f-tests  between  the  run  times  of  the 
various  matrix  exponentiation  methods  are  presented  in  Table  5.6.  The  fastest 
method  was  EE,  with  the  slowest  Uniformization,  even  though  all  three  were  fairly 
close  with  regard  to  p.  There  was  no  statistical  difference  between  the  run  times  of 
the  Uniformization  and  PSS  method,  while  the  difference  between  all  other  methods 
was  statistically  significant.  Again,  as  in  the  2-state  case,  all  three  of  the  remaining 
methods  produced  very  similar  results  for  r*,  but  the  Uniformization  method  pro¬ 
duced  a  slightly  different  value  for  A(t*).  This  suggests  the  Uniformization  method 
may  be  producing  some  slightly  inaccurate  solutions.  Based  upon  this  analysis,  the 
EE  method  is  the  most  appropriate  method  for  the  five  state  case. 

Results  from  the  truncation  experiment  are  listed  in  5.7.  The  increased  state 
dimension  and  larger  value  of  A  resulted  in  more  terms  being  truncated,  with  an 
associated  significant  decrease  in  run  time.  The  truncation  method  provided  a  17% 
improvement  over  the  baseline  method.  From  conducting  a  f-test,  this  difference  was 
highly  statistically  significant  with  a  p- value  on  the  order  of  1CU59.  For  the  baseline, 
the  mean  of  p  was  found  to  be  859.62  sec  in  the  confidence  interval  from  857.07  to 
862.18  sec.  For  the  truncation  method  the  mean  of  p  was  found  to  be  715.67  sec  in 
the  confidence  interval  from  713.81  to  717.53  sec. 
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Table  5.7  5-state  truncation  method  and  baseline  comparison. 


Baseline  Method 

Truncation  Method 

T* 

A  (r*) 

P 

°p 

5.796142578 

0.755800236 

859.63 

7.15 

5.796142578 

0.755800236 

715.67 

5.19 

5.4  7-State  Case 

The  seven-state  case  models  an  environment  with  a  moderate  number  of  states. 
The  damage  threshold  is  x  =  50  units  and  the  Poisson  shock  arrival  rate  is  A  =  1. 
The  budget  was  set  to  3.5  units  with  the  same  costs  present  in  the  other  cases.  The 
generator  matrix  is, 

-2.85  0.55  0.88  0.14  0.47  0.43  0.38 

0.73  -1.99  0.17  0.01  0.06  0.23  0.78 

0.31  0.69  -4.15  0.89  0.99  0.58  0.68 

Q  =  0.84  0.62  0.27  -3.54  0.58  0.76  0.46 

0.57  0.79  0.25  0.30  -3.01  0.53  0.57 

0.37  0.96  0.88  0.66  0.52  -4.17  0.79 

0.70  0.52  0.74  0.28  0.33  0.21  -2.79 

with  wear  rates  increasing  linearly, 


2.5  3.5  4.5  5.5  6.5  7.5  8.5 


With  a  minimum  wear  rate  of  2.5  and  x  =  50,  A  =  20  is  moderately  sized.  The  shock 
magnitude  distribution  is  uniform  over  the  interval  [0,10]  with  the  corresponding 


Laplace-Stieltjes  transform, 
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A  plot  of  the  approximate  limiting  average  availability  and  cost  over  time  for  this 
case  is  presented  in  Figure  5.3. 


Figure  5.3  Availability  (-)  and  cost  (-  -)  for  the  7-state  case. 


The  results  from  the  matrix  exponentiation  experiment  are  displayed  in  Ta¬ 
ble  5.8.  Most  notably,  the  EE  method  was  completely  unstable  and  consistently 
caused  an  exception  in  MATLAB  for  all  thirty  runs.  Even  though  Q  has  a  condition 
number  much  smaller  (on  the  order  of  105)  than  the  first  two  cases,  it  has  two  imagi¬ 
nary  eigenvalues  which  cause  the  errors.  While  the  Uniformization  and  PSS  methods 
produced  the  same  r*  to  machine  precision,  the  Taylor  series  method  produced  a 
slightly  different  value.  Furthermore,  NOMADm  performed  one  less  iteration  with 
the  Taylor  series  method  implemented.  Uniformization  and  PSS  had  fairly  close,  but 
still  significantly  different,  run  times  with  Uniformization  as  the  fastest  method.  The 
Taylor  series  method  was  significantly  slower.  The  p-values  for  these  comparisons 
are  shown  in  Table  5.9.  Based  on  this  analysis,  Uniformization  is  the  most  suitable 
matrix  exponentiation  method  for  this  case. 
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Table  5.8  7-state  matrix  exponentiation  results. 


Uniformization 

EE 

Taylor 

PSS 

T* 

1.864135742 

N/A 

1.866210938 

1.864135742 

Mr*) 

0.261806575 

N/A 

0.261735112 

0.26180678 

P 

958.1941171 

N/A 

1028.039258 

964.0604266 

oP 

6.407937487 

N/A 

4.518813375 

8.108105053 

Table  5.9  7-state  p- values  for  matrix  exponentiation  performance  tests. 


Taylor  PSS 

Uniformization 

Taylor 

3.5E-45  0.002969106 
5.9E-36 

The  truncation  method  showed  a  52%  improvement  as  presented  in  Table  5.10. 
This  difference  was  clearly  statistically  significant,  with  an  associated  p- value  from 
the  comparative  f-test  on  the  order  of  10~4'.  The  mean  value  of  p  for  the  baseline 
was  966.91  seconds  in  the  confidence  interval  of  [961.40,972.43]  seconds  and  the  mean 
p- value  for  the  truncation  method  was  465.15  seconds  in  the  interval  [464.27,466.03]. 

This  improvement  is  due  to  the  relatively  large  A  and  increased  state  size.  As 
the  state  size  grows,  the  number  of  computations  of  equation  (3.7)  grows  with  the 
square  of  the  increase,  and  any  gain  made  by  the  truncation  method  is  magnified 
with  each  computation. 


5.5  10-State  Case 

The  10-state  case  serves  as  a  more  computationally  intensive  challenge  with 
which  to  test  the  various  computational  improvements.  The  shock  arrival  rate  is 

Table  5.10  7-state  truncation  method  and  baseline  comparison. 


Baseline  Method 

Truncation  Method 

T* 

MM 

P 

Op 

1.864135742 

0.26180678 

966.9130686 

15.42015957 

1.864135742 

0.261806778 

465.1510443 

2.459362772 
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A  =  1,  and  the  critical  damage  threshold  is  x  =  60.  The  total  budget  was  set  to  2 
units.  The  following  generator  of  dimension  10  was  constructed, 


-3.86 

0.06 

0.49 

0.51 

0.07 

0.49 

0.41 

0.93 

0.37 

0.53 

0.39 

-3.37 

0.04 

0.45 

0.38 

0.50 

0.56 

0.26 

0.25 

0.55 

0.25 

0.46 

-3.93 

0.33 

0.37 

0.84 

0.27 

0.20 

0.92 

0.28 

0.35 

0.86 

0.33 

-4.67 

0.48 

0.81 

0.78 

0.05 

0.63 

0.37 

0.74 

0.86 

0.90 

0.89 

-6.18 

0.86 

0.39 

0.61 

0.88 

0.06 

0.65 

0.47 

0.31 

0.76 

0.34 

-4.30 

0.03 

0.55 

0.64 

0.54 

0.94 

0.79 

0.25 

0.88 

0.25 

0.57 

-5.41 

0.10 

0.80 

0.84 

0.83 

0.66 

0.43 

0.46 

0.58 

0.61 

0.56 

-4.71 

0.44 

0.15 

0.47 

0.00 

0.84 

0.80 

0.52 

0.10 

0.20 

0.44 

-3.55 

0.17 

0.63 

0.13 

0.18 

0.13 

0.16 

0.16 

0.09 

0.07 

0.10 

-1.65 

The  linear  wear  rates  were  designed  to  represent  environmental  states  whose 
damage  rates  are  increasing  in  a  logarithmic  fashion.  The  vector  of  rates  is 


r  = 


0.400  0.903  1.431  1.806  2.097  2.335  2.535  2.709  2.863  3.000  , 


with  the  minimum  linear  wear  rate  of  0.400  and  x  =  60,  the  maximum  lifetime  is 
A  =  150. 


The  shock  magnitude  is  Erlang  distributed  with  a  Laplace-Stieltjes  transform 


of 


Fy(u)  = 


0.5 


0.5  +  u/ 

A  plot  of  the  approximate  limiting  average  availability  and  cost  over  various  values 
of  r  is  presented  in  Figure  5.4. 


The  results  of  the  matrix  exponentiation  experiment  for  the  10-state  were  re¬ 
vealing.  As  in  the  7-state  case,  some  of  the  eigenvalues  of  Q  are  imaginary  which 
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Figure  5.4  Availability  (-)  and  cost  (-  -)  for  the  10-state  case. 


cause  the  EE  method  to  fail.  The  Uniformization  and  Taylor  series  method  both  did 
not  provide  solutions  in  the  10-state  case,  demonstrating  instability  due  to  cumu¬ 
lative  round  off  error  and  machine  tolerance  limits  as  discussed  in  chapter  4.  Only 
the  PSS  method  converged  with  a  r*  of  0.75183  and  an  associated  A(r*)  of  0.94437. 
The  average  run  time  was  31574.84  sec  with  standard  deviation  of  166.32  sec. 

The  results  of  comparing  the  truncation  to  the  baseline  for  the  10-state  case 
are  presented  in  Table  5.11.  The  truncation  method  requires  67%  less  computation 
than  the  baseline  and  is  significant  with  a  p- value  on  the  order  of  10~71.  The  mean 
value  of  p  for  the  10-state  case  is  31,764.53  seconds,  with  an  associate  confidence 
interval  of  [31,675.62,31853.43].  The  p  for  the  truncation  method  has  a  mean  of 
10,637.64  seconds  in  the  confidence  interval  of  [10602.81,10672.47]. 
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Table  5.11  10-state  truncation  and  baseline  comparison. 


Baseline  Method 

Truncation  Method 

T* 

A(C) 

P 

P * 

0.751831055 

0.944372809 

31764.53 

248.44 

0.751831055 

0.944372809 

10637.64 

97.34 

5.6  20-State  Case 


The  twenty-state  case  was  developed  to  stress  the  algorithmic  improvements 
presented  in  this  thesis  and  model  the  state-size  of  a  comprehensive,  realistic,  sce¬ 
nario.  The  critical  damage  threshold  is  x  =  100  and  the  total  maintenance  budget 
is  1.2  units.  The  20-dimension  generator  for  this  case  is  displayed  in  Tables  5.12  and 
5.13. 


Table  5.12  Twenty-state  generator  matrix  (columns  1  through  10). 


"  -10.07 

0.07 

0.60 

0.30 

0.56 

0.02 

0.81 

0.37 

0.52 

0.95  ... 

0.72 

-9.37 

0.34 

0.46 

0.70 

0.40 

0.70 

0.72 

0.14 

0.62  ... 

0.84 

0.45 

-10.38 

0.50 

0.29 

0.80 

0.44 

0.29 

0.32 

0.46  ... 

0.82 

0.21 

0.12 

-10.01 

0.99 

0.16 

0.47 

0.79 

0.86 

0.65  ... 

0.54 

0.74 

0.52 

0.94 

-10.04 

0.31 

0.88 

0.21 

0.43 

0.98  ... 

0.36 

0.58 

1.00 

0.65 

0.50 

-9.76 

0.65 

0.43 

0.06 

0.42  ... 

0.50 

0.79 

0.31 

0.95 

0.43 

0.59 

-10.01 

0.63 

0.89 

0.98  ... 

0.33 

0.08 

0.51 

0.46 

0.52 

0.05 

0.21 

-8.55 

0.26 

0.52  ... 

0.23 

0.38 

0.82 

0.27 

0.05 

0.87 

0.13 

0.25 

-8.58 

0.54  ... 

0.60 

0.20 

0.30 

0.07 

0.96 

0.73 

0.54 

0.66 

0.22 

-8.86  ... 

0.05 

0.65 

0.71 

0.19 

0.07 

0.02 

0.84 

0.13 

0.59 

0.85  ... 

0.74 

0.12 

0.87 

0.35 

0.51 

0.74 

0.24 

0.51 

0.56 

0.55  ... 

0.49 

0.51 

0.66 

0.91 

0.71 

0.29 

0.02 

0.27 

0.87 

0.82  ... 

0.99 

0.11 

0.86 

0.01 

0.62 

0.76 

0.15 

0.59 

0.22 

0.49  ... 

0.98 

0.03 

0.93 

0.99 

0.23 

0.10 

0.48 

0.27 

0.30 

0.12  ... 

0.14 

0.05 

0.94 

0.17 

0.63 

0.23 

0.88 

0.78 

0.59 

0.14  ... 

0.93 

0.18 

0.67 

0.40 

0.15 

0.49 

0.86 

0.15 

0.90 

0.78  ... 

0.02 

0.51 

0.60 

0.40 

0.31 

0.87 

1.00 

0.77 

0.52 

0.27  ... 

0.73 

0.65 

0.40 

0.79 

0.65 

0.58 

0.76 

0.52 

0.67 

0.56  ... 

0.73 

0.51 

0.13 

0.51 

0.84 

0.42 

0.89 

0.59 

0.68 

0.14  ... 
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Table  5.13  Twenty-state  generator  matrix  (columns  11  through  20). 


0.42 

0.68 

0.69 

0.56 

0.74 

0.79 

0.34 

0.51 

0.83 

0.30  ' 

0.06 

0.86 

0.97 

0.69 

0.49 

0.33 

0.02 

0.05 

0.48 

0.63 

0.80 

0.10 

0.28 

0.76 

0.65 

0.78 

0.77 

0.99 

0.34 

0.53 

0.56 

0.39 

0.17 

0.82 

0.04 

0.78 

0.46 

0.31 

0.54 

0.87 

0.77 

0.38 

0.92 

0.53 

0.33 

0.10 

0.04 

0.26 

0.56 

0.58 

1.00 

0.52 

0.26 

0.34 

0.43 

0.33 

0.10 

0.61 

0.93 

0.60 

0.09 

0.89 

0.38 

0.07 

0.10 

0.26 

0.96 

0.45 

0.27 

0.49 

0.72 

0.64 

0.64 

0.26 

0.68 

0.39 

0.61 

0.42 

0.30 

0.94 

1.00 

0.20 

0.32 

0.99 

0.03 

0.90 

0.25 

0.41 

0.61 

0.35 

0.11 

0.07 

0.19 

0.39 

0.32 

0.77 

0.93 

0.71 

0.89 

0.20 

. . .  -7.55 

0.80 

0.23 

0.79 

0.05 

0.24 

0.43 

0.10 

0.71 

0.10 

0.13 

-8.53 

0.22 

0.87 

0.37 

0.45 

0.05 

0.38 

0.20 

0.65 

0.70 

0.89 

-10.25 

0.09 

0.38 

0.27 

0.90 

0.75 

0.11 

0.61 

0.31 

0.99 

0.12 

-9.65 

0.34 

0.42 

0.98 

0.15 

0.74 

0.81 

0.06 

0.40 

0.54 

0.84 

-9.50 

0.70 

0.66 

0.66 

0.76 

0.46 

0.17 

0.34 

0.69 

0.47 

0.24 

-7.80 

0.73 

0.24 

0.28 

0.09 

0.19 

0.32 

0.72 

0.60 

0.08 

0.83 

-9.93 

0.73 

0.85 

0.10 

0.67 

0.36 

0.94 

0.21 

0.30 

0.40 

0.69 

-10.34 

0.62 

0.87 

0.39 

0.62 

0.01 

0.52 

0.18 

0.65 

0.36 

0.47 

-9.55 

0.06 

0.07 

0.60 

0.58 

0.16 

0.72 

0.15 

0.40 

0.86 

0.95 

-9.91 

The  rates  of  linear  wear  were  chosen  to  increase  quadratically  with  the  di¬ 
mension  of  Q,  as  often  occurs  when  environmental  effects  compound  to  create  dam¬ 
age.  The  following  linear  wear  rates  (to  four  decimal  positions)  were  used:  0.6444, 

1.9778,  4.2000,  7.3111,  11.3111,  16.2000,  21.9778,  28.6444,  36.2000,  44.6444,  53.9778, 
64.2000,  75.3111,  87.3111,  100.2000,  113.9778,  128.6444,  144.2000,  160.6444,  and 

177.9778.  Since  the  minimum  linear  wear  rate  is  17/18  and  x  is  100,  A  is  1800/17. 


Shocks  are  gamma  distributed  with  an  arrival  rate  of  A  =  1  and  a  Laplace- 
Stieltjes  transform, 

Fy(n)  =  (1  +  2  u)\ 


A  plot  of  the  limiting  average  availability  and  cost  over  time  for  this  case  is  presented 
in  Figure  5.5. 
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Figure  5.5  Availability  (-)  and  cost  (-  -)  for  the  20-state  case. 


The  results  from  the  20-state  matrix  exponential  experiment  show  the  com¬ 
putational  burden  incurred  by  a  generator  of  dimension  20.  Again,  the  condition 
number  of  Q  is  on  the  order  of  1016;  however,  the  real  problem  with  the  EE  method 
is  the  existence  of  imaginary  eigenvalues  that  halt  execution.  The  Taylor  series 
method,  fails  because  of  excessive  round-off  error  returning  an  excessively  large  value 
for  expm(Qt)  (the  norm  of  which  is  on  the  order  of  1051).  This  large  number  propa¬ 
gates  in  the  computations  until  the  stationary  probability  distribution  is  calculated, 
which  requires  a  matrix  division  involving  P.  In  this  case,  P  is  badly  scaled  with  a 
condition  number  of  1048,  resulting  in  false  values  for  the  stationary  transition  prob¬ 
abilities.  The  Uniformization  method  fails  for  similar  reasons.  This  analysis  shows 
that  the  PSS  method  is  best  suited  for  the  20-state  case. 

The  results  from  the  comparison  experiment  between  truncation  and  the  base¬ 
line  for  the  20-state  case  are  presented  in  Table  5.14.  Overall,  there  is  a  significant 
decrease  (with  a  p-value  on  the  order  of  10~78)  in  computational  cost  by  75%.  The 
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Table  5.14  20-state  truncation  method  and  baseline  comparison. 


Baseline  Method 

Truncation  Method 

T* 

A  (r*) 

P 

6.990850031 

0.270266593 

39367.25 

128.39 

6.990850031 

0.270266593 

9859.15 

32.47 

number  of  function  calls  was  99  in  both  cases.  The  mean  value  of  p  for  the  baseline 
case  was  39367.25  seconds  in  the  confidence  interval  of  [39,321.31,39413.20]  while 
the  mean  p  value  using  the  truncation  method  was  9859.15  seconds  in  the  interval 
[9,847.53,9,870.77], 

5. 7  Summary 

The  numerical  experiments  on  the  matrix  exponential  computation  demon¬ 
strated  the  relative  effectiveness  of  four  different  implementations.  As  a  secondary 
benefit,  agreement  in  final  solutions  provided  by  the  different  matrix  exponential 
methods  provided  confidence  in  the  results  of  the  optimization.  An  analysis  of  the 
results  over  all  cases  showed  the  Taylor  series  implementation  was  not  sufficiently 
stable  and  demonstrated,  for  this  particular  problem,  that  there  is  no  advantage  in 
any  of  the  five  cases  from  implementing  the  Taylor  series.  Also,  the  possibility  of 
imaginary  eigenvalues  demonstrated  that  the  EE  method  was  also  unsuitable,  even 
though  it  performed  the  fastest  in  the  2-  and  5-state  cases.  If,  however,  the  generator 
can  be  shown  to  be  Hermitian  (i.e.  to  have  all  real  eigenvalues),  then  the  EE  method 
may  be  an  excellent  choice. 

The  two  remaining  methods,  Uniformization  and  PSS  were  considered  on  the 
basis  of  run  time  comparisons  with  the  Uniformization  method  faster  in  the  7-state 
and  the  PSS  method  faster  in  5-state  case.  There  was  no  statistical  difference  in  the 
2-state  case.  In  the  10-  and  20-state  cases,  the  Uniformization  method  demonstrated 
some  instability.  Therefore,  based  on  the  analysis  above  and  the  recommendations 
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of  [92],  the  PSS  method  is  to  be  concluded  the  best  option  for  this  optimization 
problem.  This  was  also  the  conclusion  of  the  MathWorks  engineering  staff  in  their 
construction  of  the  expm  built-in  function. 

The  numerical  experiments  contrasting  the  truncation  method  to  the  baseline 
show  that  the  truncation  method  is  very  effective  with  savings  proportional  to  A,  the 
dimension  of  the  generator,  and  the  budget.  However,  each  characterizing  parameter 
has  an  effect  on  the  run  time.  There  is  a  general  trend,  in  that  the  truncation  method 
performs  better  as  its  associated  gains  can  be  magnified  through  increased  state  size 
and  increased  maximum  system  lifetimes.  This  matches  the  complexity  analysis  in 
that  the  algorithmic  complexity  is  0(l2 7)  with  7  directly  proportional  to  A  for  a 
fixed  r. 

Overall,  both  experiments  demonstrated  an  effective  method  for  maximizing 
limiting  average  availability  through  selecting  an  appropriate  inter-inspection  inter¬ 
val  t.  The  truncation  method  enables  the  approximate  optimal  r*  to  be  computed 
within  a  reasonable  amount  of  time.  However,  run  times  are  still  considerable  when 
maximum  lifetimes  (A)  and  the  dimension  of  Q  are  large.  Research  into  additional 
computational  techniques  could  further  reduce  this  run  time,  but  the  development  of 
a  purely  analytical  technique  could  provide  the  true  optimal  r*  with  almost  no  com¬ 
putational  burden.  This  and  other  possible  extensions  to  this  research  are  discussed 
in  the  next  chapter. 
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6.  Conclusions  and  Future  Research 


The  goal  of  this  research  was  to  present  a  methodology  to  determine  the  ap¬ 
proximate  optimum  inter-inspection  duration  (r)  that  maximizes  the  limiting  aver¬ 
age  availability  of  a  system  while  keeping  the  downtime,  inspection  and  replacement 
costs  within  an  arbitrary  budget  for  a  system  with  deterioration  due  to  environment- 
dependent  wear  and  random  shocks.  In  order  to  be  useful,  this  methodology  must 
provide  consistent  and  accurate  results  in  a  reasonable  amount  of  time.  To  this  end, 
a  secondary  goal  of  the  research  was  to  improve  the  computational  performance  of 
the  chosen  methodology. 

An  optimization  methodology  using  generalized  pattern  search  was  chosen 
and  implemented  to  compute  the  availability  measure  A(t).  Moreover,  a  numer¬ 
ical  method,  the  truncation  method ,  was  developed  which  provided  computational 
savings  most  directly  proportional  to  the  dimension  of  Q  and  the  component’s  max¬ 
imum  lifetime.  The  budget  and  other  model  parameters  also  impact  the  computa¬ 
tional  savings  of  the  truncation  method.  Several  standard  methods  were  employed 
to  improve  computational  improvement.  These  included:  vectorization,  preprocess¬ 
ing,  and  choosing  the  correct  MATLAB  operators.  In  order  to  properly  implement 
the  third  method,  a  study  of  the  matrix  exponential  computation  was  conducted  to 
determine  the  most  proper  method  for  this  research.  The  results  of  this  implemen¬ 
tation  were  presented  and  discussed  using  five  different  cases,  each  with  a  distinct 
set  of  parameters.  The  parameters  were  chosen  in  order  to  illustrate  the  means  by 
which  to  compute  the  approximate  optimal  inter-inspection  duration  and  to  evaluate 
the  associated  performance  of  the  optimization  methodology  over  a  wide-range  of 
potential  problems.  These  numerical  experiments  demonstrated  that  the  truncation 
method  presents  significant  potential  for  performance  gains  over  the  full-sum  case 
with  no  loss  of  solution  quality  (within  the  tolerance  machine  precision). 
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It  is  recommended  that  the  NOMADm  pattern  search  code  be  used,  along 
with  the  strategies  discussed  to  enhance  computational  performance  of  the  objective 
function  including  the  truncation  method.  NOMADm  was  observed  to  be  a  capable 
implementation  of  generalized  pattern  search.  Moreover,  NOMADm  is  distributed 
without  cost  under  the  GNU  General  Public  License  and  is  available  direclty  from  the 
author.  Other  recommendations  include  using  the  specific  one-dimensional  Laplace 
transform  inversion  algorithm  provided  by  Abate  and  Whitt  [5]  and  the  matrix 
exponential  should  be  computed  via  the  Pade  approximation  using  repeated  squaring 
and  scaling. 

Although  this  research  presents  a  complete  methodology  considering  degrada¬ 
tion  due  to  environment-dependent  wear  and  random  shocks,  much  work  remains. 
For  instance,  to  compute  an  approximate  value  for  r*,  the  various  parameters  de¬ 
scribing  the  environment  (shock  magnitudes  and  arrival  rate,  environment  transition 
rates,  etc.)  must  be  known.  When  these  parameters  are  unknown,  statistical  esti¬ 
mation  procedures  must  be  employed.  A  sensitivity  analysis  could  reveal  the  most 
critical  parameters  and  provide  insight  into  the  impact  that  various  environmental 
parameters  have  on  the  limiting  average  availability  and  reliability  measures.  Fur¬ 
ther  techniques  available  to  operations  research  analysts  could  determine  the  optimal 
allocation  of  resources  to  impact  the  environmental  parameters  to  maximize  long  run 
availability  or  minimize  cost. 

Additional  computational  improvement  techniques  remain  to  be  tested  for  com¬ 
puting  A(t).  In  particular,  since  the  primary  computational  expense  when  comput¬ 
ing  A(t)  is  the  population  of  the  transition  probability  matrix,  the  application  of 
current  research  [124]  in  compressing  transition  probability  matrices  could  prove 
useful.  Moreover,  while  all  computations  in  this  thesis  were  performed  in  MATLAB, 
implementing  the  methodology  in  a  compiled  language  such  as  C++  or  FORTRAN  would 
undoubtedly  improve  the  overall  run  time  by  a  significant  margin.  Further  research 
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comparing  the  results  of  this  thesis  to  meta-heuristic  methods  could  yield  computa¬ 
tionally  fast  solutions  that  also  provide  approximate  optimal  solutions. 

Additional  work  on  the  analytic  model  could  prove  fruitful.  If  a  closed-form 
expression  can  be  found  for  A(t)  in  the  time  domain,  the  true  optimal  point  could  be 
found  in  seconds.  A  pure  analytical  solution  would  also  provide  considerable  insight 
into  the  system  dynamics.  Also,  if  a  mathematically  rigorous  solution  methodol¬ 
ogy  can  be  developed  to  optimize  the  objective  function  in  the  complex  domain, 
numerical  inversion  would  not  be  necessary,  and  the  solution  could  consequently  be 
computed  much  faster  and  with  greater  precision.  Additional  insight  on  the  dynam¬ 
ics  of  the  stochastic  model  could  yield  a  maintenance  policy  more  effective  than  the 
deterministic  inspect-and-replace  policy  presented  in  this  thesis. 

There  is  also  great  potential  for  applying  this  methodology  to  an  actual  opti¬ 
mal  maintenance  setting.  Much  could  be  learned  regarding  the  applicability  of  the 
model  from  determining  environmental  transition  rates,  fitting  distributions  and  pa¬ 
rameters  and  comparing  empirical  experiments,  simulation  and  calculated  reliability 
and  availability  measures.  Moreover,  since  the  most  restrictive  assumption  is  that 
of  a  Markovian  environment,  there  is  a  great  deal  of  value  in  using  the  method  of 
Cox  [48],  who  used  phase-type  distributions  to  model  to  approximate  semi-Markov 
environment  as  a  Markovian  one. 
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Appendix  A.  Availability  Calculation 

function  [fx,cx]  =  availabilityCalc (tau) ; 

°/0availabilityCalc :  function  to  calculate  the  availability  of  a  Markovian  shock 
7,  and  wear  model  as  described  in  Kharoufeh,  et  al.  (2006). 

7. 

7,  Syntax: 

7,  fx  =  availabilityCalc  (tau) ; 

7,  [fx,cx]  =  availabilityCalc  (tau) ; 

7, 

7,  Description: 

7,  This  Matlab  program  is  used  to  compute  conditional  distributions  and 

7,  conditional  expectations  of  random  quantities  from  the  Markovian 

7,  shock  and  wear  model  in  Kharoufeh,  et  al.  (2006). 

7. 

7*  References:  Kharoufeh,  J.,  Finkelstein,  D.  ,  and  D.  Mixon  (2006) 

7,  Availability  of  inspected  systems  subject  to  Markovian  wear 

7.  and  shocks.  Technical  report.  Department  of  Operational 

7,  Sciences.  Air  Force  Institute  of  Technology. 

7. 

7*  Calls:  invt,  unif ,  expm,  ceil,  min,  ones,  zeros,  sum,  CalcLoopCnt, 

7,  getappdata,  num2str,  eye,  disp,  floor 

7. 

y^^e3|e^c^e3|e^c3)e3|e^c3|c3|e^c3|e^e^o|e^c^e3|e^c3|e^e^c3|e3|e^c3|c^e^e3|e^c^o|e^e3f;3|e^c3|e^e^c3|e^e^e3|c^c^e^e^c^;^e^c3|e3|e^e^c^e^e3|e^c^e^e^c3f;3|e^c3|c^e^e3|e^e^e3|e^c3f;^e^c^e^e^c 

7,  Copyright  (c)  2006  by  Timothy  B.  Booher 

7. - 

7,  Originally  created,  2005. 

7,  Last  modified,  1  March  2006 
7. 

7,  Author  information: 

7,  Timothy  B.  Booher,  Capt,  USAF 
7,  Air  Force  Institute  of  Technology 
7,  Department  of  Operational  Sciences 
7,  2950  Hobson  Way 

7,  Wright-Patterson  AFB,  OH  45433 
7*  Timothy.Booher@afit.edu 

7, 

7. - 

7,  VARIABLES  (only  for  the  availabilityCalc  function)  : 

7 * I N P UT * V A R I A B L E S ************************************************************** * 

7,  tau  =  the  inter-inspection  interval  and  decision  variable 

7*FUNCTI ON* VARIABLES*** ************************************************* ******** 
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42  l  1  =  overall  dimension  (number  of  states  in  Markovian  environment) 

43  /  X  =  marginal  PMF  of  R_1 

44  %  Q  =  infinitesimal  generator  matrix  for  random  environment 

45  %  r  =  vector  of  wear  rates  for  the  Markovian  environment 

46  7,  lambda  =  shock  arrival  rate 

47  7,  lim  =  the  failure  threshold  (x  in  Kharoufeh  et  al.) 

48  7*  ER_1  =  the  expected  value  (unconditional)  of  the  time  to  the  first 

49  7o  replacement  epoch 

50  7«  EF_1  =  the  expected  value  (unconditional)  of  the  first  failure 

51  7o  epoch 

52  7o  Pt  =  transition  probability  function 

53  7*  P  =  transition  probability  matrix  (of  the  embedded  DTMC) 

54  7«  p  =  stationary  distribution,  p,  of  P 

55  7*  z0  =  e_i,  vector  of  ones  with  an  1  in  the  ith  position 

56  7o  M  =  ratio  of  T_max  to  tau  (must  be  less  than  one)  the  ceiling 

57  7o  gamma  =  the  number  of  inspection  intervals  necessary  to  cover  x/r_l 

58  7*  Param.Q  =  infinitesimal  generator  matrix 

59  7«  Param.r  =  vector  of  wear  rates 

60  7o  Param.f  =  LST  of  the  c.d.f.  for  the  shock-damage  magnitudes 

61  7o  Param. lambda  =  shock  arrival  rate 

62  7o  Param. lim  =  critical  damage  threshold 

63  7o  Param. budget  =  fixed  budget  for  long-run  cost 

64  7«  Param.  cD  =  cost  of  downtime 

65  7o  Param. cl  =  cost  of  each  inspection 

66  7*  Param.  cR  =  cost  of  replacement 

67  7o  Param. L  =  x/r_l 

68  7*  Param.  states  =  dimension  of  Param.Q 

69  7«  Param.  I  =  identity  matrix  of  dimension  Param.  states 

70  7«  Param.  el  =  vector  of  elements  (ith  column  of  Param.  I) 

71  7o  P_hat  =  temp  matrix  used  to  calculate  stationary  probabilities 

72  7«  {i,k,n}  =  index  variables 

73  7.  S  =  temp  variable  to  store  a  summed  sequence 

74  3|e3|e3(c3(c3|e3(c3(e3fe3(e3|e3)e3fc3|c3fe3(c3|e3(c3|c9|e3)e3(e3te3fe3|e3fe3fe3|c3(e3|e3|e3)e3(c9|e3(e3(e3fe3(e3|e3te3(e3|e3fe3fe9|e3(e3(e3fe3fe3(c3|e3fe3|e3fe3(c3|c3fe3(e3|e3(e3(c3te3fe3(e3fe3fe3|e3fe3(e3|e3fe3(c3|c3(e 

75  7o  Load  Parameter  Information 

76  7. - 

77  Param  =  getappdata(0, ’PARAM5 ) ; 

78  disp( [’Function  call  with  \tau  =’  num2str (tau)] ) ; 

79  X  =  []  ;  P  =  []  ; 

80  7o - 

81  7o  The  code  below  is  used  to  obtain  the  transition  probability  matrix  P. 

82  1  =  Param. states ;  %  used  to  keep  things  simple 

83  gamma  =  ceil (Param. lim/ (min (Param. r) *tau) ) ; 

84 
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86 
87 


90 

91 
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94 

95 

96 

97 

98 

99 
100 
101 
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103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 
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7,  this  populates  the  transition  probability  matrix  for  the  process  that 
7,  describes  the  state  of  the  environment  at  each  replacement 
for  i=l:l 

zO=zeros(size(Param.Q(l, :))) ; 
z0(i)  =  1; 
for  k=l : 1 

for  n  =  1: gamma 

Pt  =  expm(Param. Q*n*tau) ; 

X(n)  =  invt ( ’getG5 ,z0,n*tau) ; 

PV(n)  =  Pt (i,k) ; 
if  (Param.Method(l) .Use  ==  1) 

if  ((n  >  1)  &&  ( (X(n) -X(n-l) )  <=  Param. Method (2) .Param)) 
break; 

end 

end 

end 

P(i,k)  =  sum( (X (2 : end) — X ( 1 : (end-1) ) ) * (PV (2 : end) ’))+  PV(1)*X(1); 

end 

end 

7. - 

7,  Now  compute  the  stationary  distribution,  p,  of  P 

7. - 

P_hat  =  P  -  Param. I; 

P_hat ( : ,1)  =  1; 
p  =  Param . p/P_hat ; 

7. - 


7. - 

7,  This  routine  is  used  to  obtain  ER_1 
7. - 

ER_1  =  zeros(l,l); 
for  i=l:l 


for  n=l : (gamma-1) 

S(n)=invt ( ’getG5 , Param. I (i, :) ,n*tau) ; 

end 

ER_l(i)=  tau*(gamma  -  sum(S)); 

end 


7. - 

7,  Compute  the  limiting  average  availability 
7. - 
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fx  =  -(p*Param.EF_l’)/(p*ER_l’) ; 

129 

/. - 

130 

7,  Compute  the  final  costs 

131 

•/. - 

132 

cterm  =  Param. cD* (ER_1  -  Param.EF_l)  ... 

133 

+  Param. cl*f loor (ER_l/tau)  ... 

134 

+  Param. cR*ones ( 1, 1) ;s 

135 

cx  =  (p*ctermJ )/(p*ER_l ’ )  -  Param . budget ; 

136 

return; 

137 

•/. - 

138 

i  unif ormization 

139 

y. - 

140 

function  out  =  unif (tau, lambda, P ,n,numSt at es) 

141 

pi  =  0; 

142 

ii 

143 

for  k  =  l:n 

144 

y  =  y*p*(iambda*tau)/k; 

145 

pi  =  pi  +  y; 

146 

end 

147 

out  =  exp(-lambda*tau) *pi ; 

148 

return; 

149 

150 

function  [intLoopCnt,  lambda,  P]  =  CalcLoopCnt (tau,  Q,  epsilon,  M) 

151 

K  =  0;  zi  =  1;  sigma  =  1; 

152 

lambda  =  max(-diag(Q) ) ; 

153 

eta  =  (1-epsilon) / (exp (-lambda*tau) ) ; 

154 

while  (sigma  <=  eta) 

155 

K  =  K  +  1; 

156 

zi  =  zi*(lambda*tau)/K; 

157 

sigma  =  sigma  +  zi; 

158 

end 

159 

P  =  (lambda" (-1) ) * (lambda*eye (size (Q) ) +Q) ; 

160 

intLoopCnt  =  K; 

161 

return; 

162 

163 

function  E  =  expmdemo3(A) 

164 

70EXPMDEM03  Matrix  exponential  via  eigenvalues  and  eigenvectors. 

165 

/  E  =  EXPMDEM03(A)  illustrates  one  possible  way  to  compute  the  matrix 

166 

/  exponential.  As  a  practical  numerical  method,  the  accuracy 

167 

7,  is  determined  by  the  condition  of  the  eigenvector  matrix. 

168 

•/. 

169 

"/,  See  also  EXPM,  EXPMDEM01 ,  EXPMDEM02 . 

170 
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171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 
213 


7,  Copyright  1984-2004  The  MathWorks,  Inc. 

7,  $Revision:  1. 1.6.3  $  $Date:  2004/08/16  01:37:28  $ 

[V,D]  =  eig(A) ; 

E  =  V  *  diag(exp(diag(D) ) )  /  V; 

return; 

function  E  =  expmdemo2 (A) 

7oEXPMDEM02  Matrix  exponential  via  Taylor  series. 

7,  E  =  expmdemo2(A)  illustrates  the  classic  definition  for  the 
7,  matrix  exponential.  As  a  practical  numerical  method, 

7,  this  is  often  slow  and  inaccurate. 

7. 

7.  See  also  EXPM,  EXPMDEM01 ,  EXPMDEM03 . 

7,  Copyright  1984-2003  The  MathWorks ,  Inc . 

7.  $Revision:  1. 1.6.2  $  $Date:  2004/04/10  23:24:39  $ 

E  =  zeros (size (A) ) ; 

F  =  eye(size(A)) ; 
k  =  1; 

while  norm(E+F-E, 1)  >  0 
E  =  E  +  F; 

F  =  A*F/k; 
k  =  k+1; 

end 

return; 

function  E  =  expmdemol(A) 

7oEXPMDEM01  Matrix  exponential  via  Pade  approximation. 

7.  E  =  EXPMDEMOl(A)  is  an  M-file  implementation  of  the  built-in 
7,  algorithm  used  by  MATLAB  for  the  matrix  exponential. 

7,  See  Golub  and  Van  Loan,  Matrix  Computations,  Algorithm  11.3-1. 
7. 

7.  See  also  EXPM,  EXPMDEM02 ,  EXPMDEM03 . 

7,  Copyright  1984-2003  The  MathWorks,  Inc. 

7.  $Revision:  1. 1.6.2  $  $Date:  2004/04/10  23:24:38  $ 

7,  Scale  A  by  power  of  2  so  that  its  norm  is  <  1/2  . 

[f,e]  =  log2(norm(A, ’ inf J ) ) ; 
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214  s  =  max(0,e+l) ; 

215  A  =  A/2~s; 

216 

217  1  Pade  approximation  for  exp (A) 

218  X  =  A; 

219  c  =  1/2; 

220  E  =  eye(size(A))  +  c*A; 

221  D  =  eye(size(A))  -  c*A; 

222  q  =  6; 

223  p  =  1; 

224  for  k  =  2:q 

225  c  =  c  *  (q-k+1)  /  (k*(2*q-k+l) ) ; 

226  X  =  A*X; 

227  cX  =  c*X ; 

228  E  =  E  +  cX; 

229  if  P 

230  D  =  D  +  cX; 

231  else 

232  D  =  D  -  cX; 

233  end 

234  p  =  ~p; 

235  end 

236  E  =  D\E; 

237 

238  #/0  Undo  scaling  by  repeated  squaring 

239  for  k=l:s,  E  =  E*E;  end 

240 

241  return; 


Appendix  B.  Availability  Parameters  File 


function  Param  =  availability_Param 
load  RunData; 

myCase  =  RunData. int St ate ; 

Param . Method  =  RunData. Method; 


switch  myCase 
case  {2} 

7,  2  state  case 

Param. Q  =  [-25/3  25/3; 25/3  -25/3]; 

Param. r  =  [13/12  1/4]; 

Param. f  =  inline (’ [4. 0/(4. 0+s)  4.0/(4.0+s)] ’ , ’s’) ; 

Param. lambda  =  0.50; 

Param. lim  =  1.0; 

Param. budget  =  35; 

Param. cD  =  0.5; 

Param. cl  =  1; 

Param. cR  =  5; 

case  {5} 

7,  5  state  case 

Param. Q  =  [-0.5  0.125  0.125  0.125  0.125; 


[-0.5 

0. 

125 

0. 

125 

0.125 

0. 

.  125 

0. 

.4 

-2  0. 

.4  0. 

6  0.6: 

0. 

.025 

0. 

025 

-0.1 

0.025 

0. 

.025 

0. 

.05 

0. 

05 

0. 

05 

-0.2 

0. 

.05; 

1. 

.5 

1 

1 

1. 

5  -5]  : 

Param. r 
Param. f 


=  [1  2  3  4  10] ; 

=  inline ( ’ [1  111  1] * (0 . 20/ (0 . 20+s) ) , ’ s 5 ) ; 


Param. lambda  =  0.25; 

Param. lim  =  100.0; 

Param . budget  =  0.7; 

Param. cD  =  0.5; 

Param. cl  =  1; 

Param. cR  =  5; 

case  {7} 

7.  exponential  with  7  states 


5  Param. Q 

=  [-2.8452 

0.5466 

0.8801 

0.1365 

0.4692 

0.4329 

0.3798 

6 

0.7271 

-1.9859 

0.1730 

0.0118 

0.0648 

0.2259 

0 . 7833 

7 

0.3093 

0.6946 

-4.1467 

0.8939 

0.9883 

0.5798 

0.6808 

8 

0.8385 

0.6213 

0.2714 

-3.5355 

0.5828 

0.7604 

0.4611 

9 

0.5681 

0.7948 

0.2523 

0.2987 

-3.0116 

0.5298 

0.5678 

0 

0.3704 

0.9568 

0.8757 

0.6614 

0.5155 

-4.1742 

0.7942 

1 

0.7027 

0.5226 

0.7373 

0 . 2844 

0.3340 

0.2091 

-2.7901] 

42 


Param.r  =  [1  2  3  4  5  6  7] +1 . 5 ; 

7.  uniform  0,10 

Param.f  =  inline(’[l  11111  1] * ( (exp(-s*10) -exp(-s*0) ) / (s* (10-0) ) ) ’ , ’ s ’ ) ; 


43 

44 

45  Param. lambda  =  1; 

46  Param.  lim  =  50; 

47  Param. budget  =  3.5; 

48  Param.  cD  =  0.5; 

49  Param. cl  =  1; 

50  Param. cR  =  5; 

51  case  {10} 


52  l  weibull  with  10  states 


53 

Param 

Q  =  [.•• 

54 

-3.861608355 

0.058187848 

0.494874756 

0.508179212 

55 

0.065313688 

0.486398149 

0.41364986 

0.933229859 

0.374292584 

0.527482399 

56 

0.391617223 

-3.367934425 

0.038333157 

0.452239667 

57 

0.375144658 

0.496060749 

0.560410449 

0.259379974 

0.249102724 

0.545645825 

58 

0.252783651 

0.455725962 

-3.93287784 

0.325584439 

59 

0.37352297 

0.843194083 

0.268677291 

0.204171395 

0.924875292 

0 . 284342757 

60 

0.354381927 

0.863086892 

0.327882957 

-4.669336947 

61 

0.484022397 

0.806198215 

0.784254164 

0.049208442 

0.629499293 

0.370802659 

62 

0.742977968 

0.85519697 

0.899468788 

0.886479964 

63 

-6.178942389 

0.857785596 

0.387870785 

0.60616094 

0.878308777 

0.064692602 

64 

0.650832022 

0.472255683 

0.313730488 

0.76126076 

65 

0.342061094 

-4.303955904 

0.030983621 

0.54634874 

0.641674431 

0.544809064 

66 

0.939793041 

0.78692439 

0.251676012 

0.883766335 

67 

0.252689256 

0.565730368 

-5.411183076 

0.095837436 

0.798390636 

0.836375603 

68 

0.832799125 

0.655982138 

0.432989132 

0.457406256 

69 

0.584886917 

0.611898522 

0.558558528 

-4.714868647 

0.435026039 

0.145321992 

70 

0.469977867 

3 . 98896E-05 

0.842382237 

0.799202291 

71 

0.523703594 

0.102976522 

0.200695572 

0.442948328 

-3.553446558 

0.171520257 

72 

0.629865607 

0.131237162 

0.184488992 

0.13407712 

73 

0.163419046 

0.158315917 

0.087421885 

0.066381959 

0.095957798 

-1.651165486] ; 

74  Param.r  =  [0.4  0.9031  1.4314  1.8062  2.0969  2.3345  2.5353  2.7093  2.8627  3.0000]; 

75  stri  =  ones(l,10); 

76  7,  mean  =  2 

77  7.  var  =  1 

78  7o  normally  distributed 

79  Param.f  =  inline ([’[’  int2str (stri)  J I’KO.S/CO.S+s))^’]  ,  ’s’)  ; 

80  Param. lambda  =  1; 

81  Param. lim  =  60; 

82  Param. budget  =  2; 

83  Param.  cD  =  0.5; 

84  Param.  cl  =  1; 
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85  Param.cR  =  5; 

86  case  {20} 

87  Param. Q  =  [ . . . 


88 

-10 

068 

0 

068 

0 

599 

0 

302 

0 

562 

0 

018 

0 

813 

0 

369 

0 

519 

0 

953 

0 

425 

0 

681  ... 

89 

0 

685 

0 

556 

0 

744 

0 

790 

0 

340 

0 

514 

0 

833 

0 

295; 

90 

0 

722 

-9 

370 

0 

336 

0 

460 

0 

701 

0 

404 

0 

695 

0 

716 

0 

141 

0 

624 

0 

058 

0 

856  ... 

91 

0 

974 

0 

689 

0 

489 

0 

331 

0 

019 

0 

046 

0 

476 

0 

632; 

92 

0 

841 

0 

453 

-10 

381 

0 

503 

0 

285 

0 

800 

0 

441 

0 

286 

0 

316 

0 

458 

0 

797 

0 

097  ... 

93 

0 

280 

0 

761 

0 

650 

0 

782 

0 

768 

0 

991 

0 

340 

0 

531; 

94 

0 

821 

0 

211 

0 

118 

-10 

014 

0 

987 

0 

160 

0 

471 

0 

789 

0 

861 

0 

648 

0 

560 

0 

392  ... 

95 

0 

174 

0 

824 

0 

042 

0 

777 

0 

458 

0 

310 

0 

542 

0 

871; 

96 

0 

544 

0 

735 

0 

520 

0 

944 

-10.037 

0 

312 

0 

881 

0 

205 

0 

433 

0 

979 

0 

771 

0 

376  ... 

97 

0 

922 

0 

528 

0 

335 

0 

105 

0 

044 

0 

256 

0 

562 

0 

584; 

98 

0 

364 

0 

577 

0 

995 

0 

650 

0 

501 

-9 

756 

0 

645 

0 

430 

0 

059 

0 

421 

0 

995 

0 

516  ... 

99 

0 

260 

0 

342 

0 

433 

0 

329 

0 

100 

0 

610 

0 

925 

0 

603; 

100 

0 

501 

0 

789 

0 

313 

0 

946 

0 

430 

0 

591 

10 

015 

0 

629 

0 

886 

0 

976 

0 

090 

0 

894  ... 

101 

0 

375 

0 

068 

0 

099 

0 

261 

0 

957 

0 

447 

0 

274 

0 

488; 

102 

0 

330 

0 

083 

0 

506 

0 

464 

0 

516 

0 

054 

0 

212 

-8 

547 

0 

256 

0 

519 

0 

721 

0 

642  ... 

103 

0 

645 

0 

264 

0 

681 

0 

388 

0 

609 

0 

418 

0 

296 

0 

943; 

104 

0 

231 

0 

376 

0 

820 

0 

271 

0 

047 

0 

869 

0 

126 

0 

246 

-8 

582 

0 

543 

0 

997 

0 

203  ... 

105 

0 

316 

0 

994 

0 

027 

0 

897 

0 

247 

0 

412 

0 

613 

0 

347; 

106 

0 

605 

0 

201 

0 

299 

0 

069 

0 

965 

0 

728 

0 

536 

0 

657 

0 

224 

-8 

860 

0 

108 

0 

069  ... 

107 

0 

192 

0 

385 

0 

321 

0 

767 

0 

930 

0 

714 

0 

891 

0 

200; 

108 

0 

047 

0 

646 

0 

709 

0 

187 

0 

066 

0 

017 

0 

840 

0 

132 

0 

591 

0 

849 

-7 

546 

0 

800  ... 

109 

0 

233 

0 

786 

0 

052 

0 

244 

0 

430 

0 

103 

0 

714 

0 

098; 

110 

0 

737 

0 

124 

0 

869 

0 

349 

0 

510 

0 

744 

0 

240 

0 

513 

0 

564 

0 

546 

0 

133 

-8 

526  ... 

111 

0 

220 

0 

870 

0 

370 

0 

455 

0 

049 

0 

380 

0 

202 

0 

652; 

112 

0 

488 

0 

508 

0 

656 

0 

914 

0 

705 

0 

292 

0 

022 

0 

274 

0 

866 

0 

815 

0 

704 

0 

892  ... 

113 

-10 

250 

0 

085 

0 

384 

0 

272 

0 

905 

0 

752 

0 

108 

0 

606; 

114 

0 

986 

0 

106 

0 

862 

0 

014 

0 

625 

0 

756 

0 

148 

0 

587 

0 

220 

0 

486 

0 

306 

0 

992  ... 

115 

0 

124 

-9 

647 

0 

336 

0 

423 

0 

978 

0 

149 

0 

738 

0 

811; 

116 

0 

978 

0 

029 

0 

927 

0 

988 

0 

227 

0 

104 

0 

476 

0 

268 

0 

301 

0 

116 

0 

064 

0 

401  ... 

117 

0 

538 

0 

845 

-9 

503 

0 

699 

0 

658 

0 

661 

0 

762 

0 

462; 

118 

0 

142 

0 

051 

0 

939 

0 

174 

0 

631 

0 

227 

0 

878 

0 

784 

0 

587 

0 

143 

0 

174 

0 

341  ... 

119 

0 

689 

0 

466 

0 

238 

-7 

801 

0 

727 

0 

237 

0 

279 

0 

094; 

120 

0 

929 

0 

183 

0 

666 

0 

403 

0 

153 

0 

488 

0 

863 

0 

149 

0 

902 

0 

785 

0 

193 

0 

317  ... 

121 

0 

720 

0 

599 

0 

081 

0 

825 

-9 

931 

0 

728 

0 

845 

0 

103; 

122 

0 

016 

0 

510 

0 

598 

0 

399 

0 

311 

0 

873 

0 

996 

0 

773 

0 

524 

0 

269 

0 

670 

0 

364  ... 

123 

0 

939 

0 

208 

0 

301 

0 

403 

0 

686 

-10.340 

0 

624 

0 

875; 

124 

0 

728 

0 

648 

0 

405 

0 

786 

0 

649 

0 

579 

0 

757 

0 

520 

0 

669 

0 

557 

0 

387 

0 

616  ... 

125 

0 

008 

0 

517 

0 

180 

0 

652 

0 

363 

0 

470 

-9 

550 

0 

061; 

126 

0 

734 

0 

510 

0 

134 

0 

506 

0 

836 

0 

421 

0 

886 

0 

587 

0 

677 

0 

136 

0 

070 

0 

597  ... 

127 

0 

582 

0 

158 

0 

716 

0 

152 

0 

401 

0 

858 

0 

955 

-9 

915;] 

; 
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128 


r  =  [0.4444  1.7778  4.0000  7.1111  11.1111  16.0000  21.7778  28.4444]; 

129  r  =  [r  36.0000  44.4444  53.7778  64.0000  75.1111  87.1111  100.0000]; 

130  Param.r  =  [r  113.7778  128.4444  144.0000  160.4444  177.7778]; 

131  Param.r  =  Param.r  +  0.2; 

132  7.  gamma  distributed  shock  distribution 

133  °/0theta  =  1; 

134  °/0k  =  7 

135  stri  =  ones (1,20); 

136  Param.f  =  inline ([’[’  int2str (stri)  ’] *(l+2*s) ~4’] , ’s’) ; 

137  Param. lambda  =  1; 

138  Param. lim  =  100; 

139  Param. budget  =  1.2; 

140  Param.  cD  =  0.5; 

141  Param. cl  =  1; 

142  Param.  cR  =  5; 

143  otherwise 

144  erstr  =  ’7»g  is  not  a  specified  state.’; 

145  error (erstr ,myCase) 

146  r  =  0; 

147  end 

148  Param. L  =  Param. lim/min (Param. r) ; 

149  Param. states  =  length (Param. r) ; 

150  Param. R_D  =  diag (Param . r) ; 

151  Param. I  =  eye (Param. states) ; 

152  Param. el  =  ones(Param. states, 1) ;  %  Create  column  vector  of  ones 

153  7,  now  perform  one  set  of  calculations  for  the  invt  function 

154  m=ll;  c=[];  ga=8;  A=ga*log(10) ; 

155  Param.  mm=2~in; 

156  Param. c  =  [1  11  55  165  330  462  462  330  165  55  11  1]; 

157  Param. u  =  exp (A/2) /Param. lim; 

158  Param. x  =  A/(2*Param. lim) ; 

159  Param. h  =  pi/Param. lim; 

160  Param. su=  zeros (m+2) ; 

161  7,  calculate  the  p-values 

162  Param. p  =  zeros (1 , Param. states) ; 

163  Param. p (1 , Param. states)  =  1; 

164  7*  calculate  the  numerator 

165  EF_1  =  zeros (1 , Param. states) ; 

166  for  i  =  1 : Param. states 

167  EF_l(i)  =  invt ( ’getE’ , Param. I (i ,:), 1 . 0 ,  Param); 

168  end 

169  Param. EF_1  =  EF_1; 

170  return 
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Appendix  C.  Laplace  Inversion  Code 

1  ^C3|e3(c^e9|e3(e3(e3)e3(e3|e3)e3fe3|c3fe^c3|c3(e3|e3|e3(e3fc3te3fe3|c3fe3fe3|e3(e^e3|e3(e^c9|e3(e3(c3|e3fe3|e3te3(e3|e3(e^e3|e3(e3(c3ie3fe3(c3te3(e3|e3(e^c3|c3(e3(c3|e3(e3(c3te3fe3(e3fe3fe3|e3fe3(c3|c3fe3(c3|c3(e^e 

2  7.  invt:  accomplishes  numerical  inversion  of  laplace  transforms 

3  l - 

4  7,  Called  by:  IFS 

5  7«  Calls:  feval,  dot,  zeros 

6  7.  VARIABLES: 

7  7«  EG  =  full  array  (n  x  1) 

8  7o  zO  =  individual  item  to  trim  (scalar) 

9  7o  ntau  =  trimmed  array 

10  7*  c  =  index  of  individual  item  to  trim  (scalar) 

11  7«  u  =  length  of  the  full  item  to  trim 

12  7o  rho,qx,m,c,ga,mm,ntr,u,x,h,su,sm,f 1  =  algorithmic  parameters 

13  7o  k  =  index  variable 

14  7o  j  =  sqrt(-l) 

15  afcslesfcslesleslestealesfeslealesfeslesfesleslestesleslesfestcsteslesleafesfeslcsfesleslesfcsleslesfestestesfeslestesfeslcstesleslesfestcaiesfestcstesfeslcaleilcslcslestcslesfestcstesfestcaiesfeslestestcslesfeilcsles^sfesfe^fistcslesfe 

16  7«  References: 

17  7o  Abate,  J.  and  W.  Whitt  (1995).  Numerical  inversion  of  the 

18  7*  Laplace  transform  of  probability  distributions.  ORSA 

19  7«  Journal  on  Computing,  7,  36-43. 

20  function  fl  =  invt (EG, zO, ntau,  Param) 

21  if  (nargin  <  4) 

22  Param  =  getappdata(0 , ’PARAM’ ) ; 

23  end 

24  ntr=15; 

25  sm=  f eval (EG, Param. x,z0, ntau, Param. f (Param.x) , Param) /2; 

26  for  k=l:ntr 

27  sm  =  sm  +  ((-1) ~k)* . . . 

28  f  eval (EG , Param . x+k*Param . h*  j , zO , nt au , Param . f (Param . x+k*Param . h*  j ) , Param) ; 

29  end 

30  Param. su(l)=sm; 

31  for  k=l:12 

32  n  =  ntr+k; 

33  Param. su(k+l)  =  Param. su(k)  +  ( (-1) ~n) * . . . 

34  f  eval (EG , Param . x+n*Param . h*  j , zO , nt au , Param . f (Param . x+n*Param . h* j ) , Param) ; 

35  end 

36  fl  =  Param. u*(dot (Param. c , Param. su(l : 12, 1) O )/Param. mm; 

37  return; 

38  7o - 

39  7o  getE 

40  7o - 

41  function  out_getE  =  getE(s,zO,ntau,f_eval,  Param) 
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42  7.  R_D  =  diag(r) ; 

43  F  =  diag(f_eval) ; 

44  °/0  I  =  eye  (size  (Q) ) ; 

45  7,  el  =  ones (size (Q , 1) , 1) ;  7,  Create  column  vector  of  ones 

46  if  (nargin  <  4) 

47  Param  =  getappdata(0 , ’PARAM’ ) ; 

48  end 

49  7«  Compute  the  conditional  expectation  values  of  F_l. 

50  z  =  (l/s)*(zO*inv(s*Param.R_D-Param.Q-Param.lambda*(F-Param.I)))*Param.el*ntau; 

51 

52  7o  The  desired  value  is  the  real  part  of  z. 

53  out_getE  =  real(z); 

54  return; 

55  7. - 

56  7«  getG 

57  7. - 

58  function  out_getG  =  getG(s,zO,ntau,f_eval,  Param) 

59  F  =  diag(f_eval) ; 

60  if  (nargin  <  5) 

61  Param  =  getappdata(0 , ’PARAM5 ) ; 

62  end 

63  A2  =  expm ( (Param. Q+Param. lambda* (F-Param. I ) -s*Param. R_D) *nt au) ; 

64 

65  7«  Compute  the  cdf  values  for  given  t  (ntau)  . 

66  z  =  (l/s)*(l-zO*A2*Param.  el) ;  7«  get  the  cdf 

67 

68  7o  The  desired  value  is  the  real  part  of  z. 

69  out_getG  =  real(z); 

70  return; 
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function  [BestF,BestI  ,RunStats  ,RunSet]  =  mads_batch 
°/oMADS_BATCH  Sets  up  and  runs  the  MADS  algorithm  without  a  GUI. 

7, 

7,  Syntax: 

7,  mads_batch 

7. 

7,  Description: 

7,  This  function  serves  as  a  GUI-free  alternative  to  NOMADm  in  setting 

7,  up  an  optimization  problem,  setting  various  algorithm  parameters  and 

7,  user  options,  and  calling  the  MADS  optimizer.  It  first  sets  all  of  the 

7,  variables  to  their  default  values,  which  are  clearly  stated  in  the 

7,  MADS_DEFAULTS  file.  To  change  a  variable  from  its  default  value,  the 

7,  user  must  add  a  statement  to  this  file  to  do  so.  Some  variable 

7,  statements  are  included  here  for  convenience,  which  can  be  change 

7,  manually. 

7. 

7.  See  also  MADS_DEF AULTS ,  MADS 

y^^e3|e^c^e3|e3|c3)e3|e^c3|c^e^c3|e^e3f;3|c^c^c3|e^c^;^e^c3|o|c^c3|e^e^e3|c^c^e3|c^c3)o|e^c3|c3ie^cHe3le:le3le:lcHi:3lc:4c:f::ie:4C3le3le3lC3le3(e3(e3le:lc,f::iC3lC3f:3ic^Hoiof=He3le,f:3le:lC3f:3lc:lc,f::ic:lc 

7,  Copyright  (c)  2001-2005  by  Mark  A.  Abramson 
7. 

7,  This  file  is  part  of  the  NOMADm  software  package. 

7. 

7*  NOMADm  is  free  software;  you  can  redistribute  it  and/or  modify  it  under  the 
7,  terms  of  the  GNU  General  Public  License  as  published  by  the  Free  Software 

7,  Foundation;  either  version  2  of  the  License,  or  (at  your  option)  any  later 

7,  version. 

7. 

7,  NOMADm  is  distributed  in  the  hope  that  it  will  be  useful,  but  WITHOUT  ANY 
7,  WARRANTY;  without  even  the  implied  warranty  of  MERCHANTABILITY  or  FITNESS 

7,  FOR  A  PARTICULAR  PURPOSE.  See  the  GNU  General  Public  License  for  more 

7,  details . 

7. 

7,  You  should  have  received  a  copy  of  the  GNU  General  Public  License  along 

7,  with  NOMADm;  if  not,  write  to  the  Free  Software  Foundation,  Inc., 

7,  59  Temple  Place,  Suite  330,  Boston,  MA  02111-1307  USA 

7. - 

7,  Originally  created,  2001. 

7,  Last  modified,  1  March  2006  by  Tim  Booher 
7, 

7,  Author  information: 
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42  7.  Mark  A.  Abramson,  LtCol,  USAF,  PhD 

43  7,  Air  Force  Institute  of  Technology 

44  7*  Department  of  Mathematics  and  Statistics 

45  7«  2950  Hobson  Way 

46  7o  Wright-Patterson  AFB,  OH  45433 

47  7.  (937)  255-3636  x4524 

48  7*  Mark .  Abramson@af  it .  edu 

49  %5)t5(e3|e5|c5(ej)e5|es|c3)e5(es|sj|s5(es|c5|t5(ej)s5|e5(ej)e5|cs|ej|e5(es|c3)c5(es(e5|e5(ej)c5|c5(ej|e5|cs|cj|e5(es|sj|c5(es|e5|c5(ej)c5|es|cj)s5|cs|cj|s5(es|c5|c5(e3|c5|e5(ej|e5|cs|ej|e5|es|cj)c5(es|e5|e5(es|c5|c5(ej)e5|es|cj)e5|es|cj|e 

50 

52  7o  Calls:  mads_def aults ,  mads,  <  user  initial  points  file  > 

53  7*  Variables : 

54  7o  Defaults  =  structure  of  MADS  default  values  (see  mads_def aults) 

55  7o  Options  =  structure  for  options  settings  (see  mads_def aults) 

56  7o  problemPath  =  location  of  user  problem  files 

57  7o  Problem  =  structure  of  data  for  optimization  problem 

58  7«  newPath  =  logical  indicating  if  path  is  not  the  Mat  lab  path 

59  7o  iterateO  =  structure  of  data  for  the  initial  iterate  (see  mads) 

60  7o  BestF  =  final  best  feasible  solution  found 

61  7o  BestI  =  final  least  infeasible  solution  found 

62  7o  RunStats  =  structure  of  MADS  Run  statistics  (see  mads) 

63  7^^^^^^^^^^3^^c3ic3fc^i:3lc:lc3le3le:ic^i:^c^c3ic3fc^c3le^c:^3le:le^i:3lc:lc^i:^c^c3ie^c:lc3lc3^,ic3lc:ic^i:3fc^c^i:^c:lc3lc^c^i:3le^c^i:3lc:ic^i:^c^c^i:3fc^c3ic^c3^3i(^c^e3lc:lc^i:^c^c^i:3ie^c 

64 

65  7o  Set  Options  to  their  default  values 

66  clear  variables 

67  Defaults  =  mads_def aults ( ’Truth’ ) ; 

68  Options  =  Defaults . Options ; 

69  Problem. nameCache  =  ’CACHE’; 

70 

71  7o  Specify  Problem  Files 

72  problemPath  =  pwd; 

73  Problem. File .F  =  ’ availabilityCalc ’ ;  7*  functions  file 

74  Problem. File . 0  =  ’ availabilityCalc_Omega’ ;  %  linear  constraints  file 

75  Problem. File .X  =  ’ availabilityCalc_X’ ;  7*  closed  constraints  file 

76  Problem. File . I  =  ’ availabilityCalc_xO ’ ;  7*  initial  points  file 

77  Problem. File .N  =  ’ availabilityCalc_N’ ;  7*  discrete  neighbor  file  (MVP  only) 

78  Problem. File .P  =  ’ availabilityCalc_Param’ ;  7*  parameter  file 

79  Problem. File. C  =  ’ availabilityCalc_Cache .mat ’ ;  7*  previously  created  Cache  file 

80  Problem. File. S  =  ’ availabilityCalc_Session .mat ’ ;  7*  previously  created  Session  file 

81  Problem. File .H  =  ’availabilityCalc_History.txt’;  %  iteration  history  text  file 

82  Problem. fType  =  ’M’ ;  i  type  of  functions  file  {M,F,C} 

83  Problem. nc  =  1;  "/  number  of  nonlinear  constraints 

84 
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85  7,  Set  the  path,  and  load  any  user-provided  problem  parameters 

86  cwd  =  pwd;  cd(problemPath) ; 

87  if  (exist (Problem. File. P, ’file’)  ==  2) 

88  Problem. Param  =  feval (Problem . File . P) ; 

89  setappdata(0, ’PARAM’ , Problem. Param) ; 

90  end 

91 

92  /  Specify  Choices  for  SEARCH 


93 

Options . Search(l) .type 

=  ’None ’ ; 

7.  For  choices,  see  mads_def aults 

94 

Options . Search(l) .niter 

=  0; 

7.  Number  of  iterations  for  Search 

#i 

95 

Options . Search(l) .nPoints  =  0; 

7.  Number  of  poll  or  sample  points 

96 

Options . Search(l) .file 

=  >  >  • 

7#  filename  must  include  the  full  path 

97 

Options . Search(l) .local 

=  0; 

7#  flag  to  turn  on  trust  region 

98 

Options . Search(l) .merit 

=  0; 

7o  flag  to  penalize  clustered  data 

99 

Options . Search (2) .type 

=  ’None ’ ; 

7.  For  choices,  see  mads_def aults 

100 

Options . Search(2) .niter 

=  Inf; 

7.  Number  of  iterations  for  Search 

#2 

101 

Options . Search(2) .nPoints  =  0; 

7.  Number  of  poll  or  sample  points 

102 

Options . Search(2) .file 

=  > }  • 

7#  filename  must  include  the  full  path 

103 

Options . Search(2) .local 

=  0; 

7.  flag  to  turn  on  trust  region 

104 

Options . Search(2) .merit 

=  0; 

7o  flag  to  penalize  clustered  data 

105 

Options .nSearches 

=  2; 

106 

Options . SurOptimizer 

=  ’fmincon’ ; 

107 

108 

7,  Specify  Choices  for  POLL 

109 

Options .pollStrategy 

=  ’Standard_2n’ 

7.  For  choices,  see  mads_def aults 

110 

Options .pollOrder 

=  ’Consecutive’ 

7.  For  choices,  see  mads_def aults 

111 

Options .pollCenter 

=  0; 

7.  Poll  around  n-th  filter  point 

112 

Options .pollComplete 

=  0; 

7.  Flag  for  complete  polling 

113 

114 

7,  Specify  Termination  Criteria 

115 

Options .Term . delta 

=  le-4; 

7#  minimum  mesh  size 

116 

Options . Term . niter 

=  Inf; 

7.  maximum  number  of  iterations 

117 

Options . Term . nFunc 

=  50000; 

7.  maximum  number  of  function  evals 

118 

Options .Term. time 

=  Inf; 

7.  maximum  CPU  time 

119 

Options . Term . nFails 

=  Inf; 

7.  max  number  of  consecutive  Poll  fails 

120 

121 

7,  Choices  for  Mesh  Control 

122 

Options .deltaO 

=  i; 

7oProblem. Param. L/4;  %  initial 

mesh  size 

123 

Options . deltaMax 

=  Inf; 

7o  bound  on  how  coarse  the  mesh  can  get 

124 

Options .meshRef ine 

=  0.5; 

7.  mesh  refinement  factor 

125 

Options .meshCoarsen 

o 

T— 1 

II 

7.  mesh  coarsening  factor 

126 

127 

7,  Choices  for  Filter  management  (for  problems  with  nonlinear  constraints) 
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128  Options  .hmin  =  le-4;  l  minimum  infeasible  point  h-value 

129  Options.hmax  =1.0;  °/«  maximum  h-value  of  a  filter  point 

130 

131  7,  Choices  for  EXTENDED  POLL  (for  MVP  problems) 

132  Options . ePollTriggerF  =0.01;  i  f-value  Extended  Poll  trigger 

133  Options . ePollTriggerH  =0.01;  %  h-value  Extended  Poll  trigger 

134 

135  7,  MADS  flag  parameter  values 

136  Options . loadCache  =  1;  7#  load  pre-existing  Cache  file 

137  Options . countCache  =1;  l  count  Cache  points  as  function  calls 

138  Options  .runStochastic  =0;  7«  runs  problem  as  a  stochastic  problem 

139  Options . scale  =  2;  7«  scale  directions  using  this  log  base 

140  Options .useFilter  =0;  l  filter  (0=none,  l=multi-pt,  2=2-pt) 

141  Options .degeneracyScheme  =  ’random’;  7«  scheme  for  degenerate  constraints 

142  Options . removeRedundancy  =0;  l  discard  redundant  linear  constraints 

143  Options . computeGrad  =0;  7«  compute  gradient,  if  available 

144  Options . saveHistory  =0;  l  saves  MADS  performance  to  text  file 

145  Options  .plotHistoryl  =0;  7«  plot  MADS  performance 

146  Options .plotHistory2  =0;  l  plot  MADS  performance  real-time 

147  Options .plotFilter  =0;  )  plot  the  filter  real-time 

148  Options  .plotColor  =  ’k’ ;  7«  color  of  history  plot 

149 

150  7«  Set  up  figure  handles  for  real-time  plots 

151  if  (Problem. nc  ==  0) 

152  Options .plotFilter  =  0; 

153  end 

154  if  (Options .plotFilter) 

155  figure;  Options . fplothandle  =  gca; 

156  end 

157  if  (Options .plotHistory2) 

158  figure;  Options .hplothandle  =  gca; 

159  end 

160 

161  7o  Get  the  initial  iterates  and  call  the  optimizer 

162  iterateO  =  feval (Problem. File . I) ; 

163  [BestF,BestI ,RunStats ,RunSet]  =  mads (Problem, iterateO , Options) ; 

164 

165  7o  Perform  any  user-defined  post-processing  (must  have  argument) 

166  if  (exist (Problem. File. P, ’file’)  ==  2)  &&  (nargin (Problem. File .P)  <  0) 

167  Param  =  feval(Problem.File.P,BestF) ;  setappdata(0, ’PARAM’ ,Param) ; 

168  end 

169  cd(cwd) ; 

170  return 
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Appendix  E.  GPS  Constraints  File 
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y^^o|c^c^e3|e^c^e^e^c3|o|e^c3|c3|e^e3|e^c^e3|c^c^c3|e^c3|o|e^:3|c9|c^e3|c^c^e3|c^c3|o|e^:3|oie^:He3(e:le3le:lcHe3fe:lc3(e>fe:4cHe3le:lC3le3le3(C3le3lc:(c:ic:4C3f:3ie:lC3lo|eHe3le 

7.  availability_Omega:  User-supplied  function  for  defining  Omega,  based  on  p. 

7. - 

7,  Variables: 

7,  A  =  Coefficient  matrix  for  bound  and  linear  constraints 

7,  1  =  Lower  bounds  for  A*x  for  any  iterate  x 

7,  u  =  Upper  bounds  for  A*x  for  any  iterate  x 

^^o|c^c^e3|e^c3)e3)e^c3|e3|e^e3|c3|e^e3|e^e4;3|e^c^o|c^e3|e3|e^:3|e3|e3|e3|c^c3f;3|c^e3|o|e^:3|ole:4eHe:(eHe3le:4C3ii:3le:4C3(e9ie:4C3le,(e:le3le3lc,ii:3le3lc:(C9fc:t:3f:3le3|cHe:lcHe3le 

function  [A,l,u]  =  example_Omega(n) ; 

Par am  =  getappdata(0, ’PARAM’ ) ; 

A  =  eye (n) ; 

1  =  [0.01] ; 
u  =  [Param. L] ; 
return; 
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Appendix  F.  Initial  Iterates  File 

1  function  iterate  =  availabilityCalc_xO; 

2  l  iterate. p  is  a  vector  containing  the  values  of  the  continuous  variables 

3  /  iterate. x  is  a  cell  array  containing  the  values  of  the  categorical 

4  Param  =  getappdata(0, ’PARAM’ ) ; 

5  iterate(l).x  =  Param. L/2;  #/«  start  in  the  middle  of  the  lifetime 

6  iterate(l).p  =  {};  °/0  null 

7  return; 
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Appendix  G.  MADS  Batch  Script 

function  myoutput  =  CallMadsBatch(intSt , i ,rnum,  origDir) 

70CallMadsBatch:  function  to  set  parameters  and  call  a  mads  batch  run 

7, 

7,  Syntax: 

7,  CallMadsBatch (States ,  Method,  RunID,  OriginalDirectory)  ; 

7. 

7,  Description: 

7,  This  Mat  lab  program 
7. 

7*  Calls:  mads_batch,  int2str,  disp,  save,  tic,  toe 
7. 

3(e3|e3|e3)e3|e3|e3tc3|c3|e3|e3fc3|e3|e3fe3fc3|e3|e3(e3|e3|e3)c3fe3|c3|e3tc3(c3|e3fc3tc3|e3|e3(e3|e3|c3(c3|c3|c3|e3fc3|e3|e3(e3tc3|e3|c3(e3|c3|c3te3fc3|e3|e3fe3|e3|e3fe3(c3|e3|e3tc3{c3|e^c3fc3|e3|e3fe3tc3|e3|e^c3|e3|e^c3fc3|c^c3fc3(c 

7,  Copyright  (c)  2006  by  Timothy  B.  Booher 

7. - 

7,  Originally  created,  2005. 

7,  Last  modified,  14  January  2006 

7. 

7,  Author  information: 

7*  Timothy  B.  Booher,  Capt,  USAF 
7,  Air  Force  Institute  of  Technology 
7,  Department  of  Operational  Sciences 
7,  2950  Hobson  Way 

7,  Wright-Patterson  AFB,  OH  45433 
7*  Timothy.Booher@afit.edu 

y^^c3|c^c^e3|e^e^e^e^c3|o|e^c3|c^e3f;3|e^e3|e3|e^c3)c^e^e3|c3|e^:3|o|e^e3|c^c3f:3|c^c^o|e^:3|ofe:4eHe3(e3le3le:4e3f:3le:t:3(C9fe:4cHe:(e:le3le9fe3f:He3lC3(C3fe:4c3f:3fe3lcHo|oii:3le:(e:ie3le:lc:ii::ic:lc,f:9fe:t: 

7, 

7.  VARIABLES  (only  for  the  CallMadsBatch  function)  : 

7* INPUT* VARI ABLES ************************************************************** * 

7,  tau  =  the  inter-inspection  interval  and  decision  variable 

7.  intSt  =  the  state 

7,  i  =  the  method  parameter  case 

7.  rnum  =  a  unique  id  number  to  describe  the  current  run 

7,  origDir  =  the  directory  with  the  location  of  the  original  run  files 

%*FUNCTI ON* VARIABLES*** ************************************************* ******** 

7.  RunDat a. Method (n)  .Use  =  boolean  (1,0)  that  states  if  the  method  should 

7,  be  used  in  case  ’n’ 

7,  RunDat  a.  Method(n)  .Param  =  parameter  for  a  particular  method 

7*  BestF,BestI ,RunStats ,  RunSet  =  mads  output  variables  (see  mads.m) 

7.  runTime.time  =  my  record  of  the  program  run  time 

y^^e3|e^c^e3|e^c3f:3|e^c3|c^e^c3|e^e3f;3|e^c^e3|e^c3f;^e^c3|o|e^c3|c3|e3f;3|e^c^c3|e^c^o|e^c3|c^e^c3|c^e^e^c^c3|o|e^c^e^e^c3|e^e^e^c^e3|o|e^c3f;^e^c3f;^e^c3|c^e^e3|e^e^o|e^c 

7,  set  parameter  file 
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42  RunData. intState  =  intSt; 

43  %  Method  1  truncation 

44  %  Method  2  Matrix  exponential  via  Unif ormization 

45  7,  Method  3  Matrix  exponential  via  eigenvalues  and  eigenvectors 

46  7,  Method  4  Matrix  exponential  via  Taylor  series. 

47  7o  Method  5  Matrix  exponential  via  Pade  approximation. 

48  7o  runl  =  baseline  I  run2  =  trunc  I  run3  =  unif  I  run4  =  eig  I  run5  =  Taylor  I  run6  =  Pade 

49  if  (nargin  >  3) 

50  path (path, or igDir) 

51  end 

52  blnUseFlops  =  true; 

53  if  blnUseFlops 

54  if  (exist (’C:\Program  Files\ICL\WinPAPI\PAPI  MATLAB  Support\f lops . dll 5) ~=3) 

55  path (’ C:\Program  Files\ICL\WinPAPI\PAPI  MATLAB  Support\ 5 ,path) ; 

56  disp(5path  added  for  FLOPS5); 

57  else 

58  disp(5path  addition  not  needed  —  FLOPS  in  use5); 

59  end 

60  end 

61  7o  save  structure  in  .mat  file 

62  switch  i 

63  case  1 

64  7o  RunData  —  contains  information  on  each  run 

65  7.  *$*$*$*  RUN  1:  Baseline  *$*$*$* 

66  7#  1:  gamma  scaling  I  2:  truncation  I  3:  second  loop  tol  I  4:  unif  ormization 


67  °/0  *****  what  methods  do  we  want  to  run? 

68  RunData. Method (1) .Use  =  false; 

69  RunData. Method (2) .Use  =  false; 

70  RunData. Method (3) .Use  =  false; 

71  RunData. Method (4) .Use  =  false; 

72  RunData. Method (5) .Use  =  false; 

73  7o  *****  what  parameters  needed  for  each  method? 

74  RunData. Method (1) .Par am  =  0; 

75  RunData. Method(2) .Param  =  0; 

76  RunData. Method (3) .Param  =  0; 

77  RunData. Method (4) .Param  =  0; 

78  RunData. Method (4) .Param  =  0; 

79  case  2 


80  7.  *$*$*$*  RUN  2:  Truncation  *$*$*$* 

81  7o  *****  what  methods  do  we  want  to  run? 

82  RunData. Method (1) .Use  =  true; 

83  RunData. Method (2) .Use  =  false; 

84  RunData. Method (3) .Use  =  false; 
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85 

RunData. Method (4) .Use  =  false; 

86 

RunData. Method (5) .Use  =  false; 

87 

7,  *****  what  parameters  needed  for  each  method? 

88 

RunData. Method(l) .Par am  =  le-4; 

89 

RunData. Method(2) .Param  =  0; 

90 

RunData. Method (3) .Param  =  0; 

91 

RunData. Method (4) .Param  =  0; 

92 

RunData. Method (5) .Param  =  0; 

93 

case  3 

94 

7,  *$*$*$*  RUN  3:  Unif ormization  *$*$*$* 

95 

7,  1:  gamma  scaling  I  2:  truncation  I  3:  second  loop  tol  I  4:  unif  ormization 

96 

%  *****  what  methods  do  we  want  to  run? 

97 

RunData. Method (1) .Use  =  false; 

98 

RunData. Method (2) .Use  =  true; 

99 

RunData. Method (3) .Use  =  false; 

100 

RunData. Method (4) .Use  =  false; 

101 

RunData. Method (5) .Use  =  false; 

102 

7,  *****  what  parameters  needed  for  each  method? 

103 

RunData. Method(l) .Param  =  0; 

104 

RunData. Method (2) .Param  =  le-2; 

105 

RunData. Method (3) .Param  =  0; 

106 

RunData. Method (4) .Param  =  0; 

107 

RunData. Method (5) .Param  =  0; 

108 

case  4 

109 

7,  *$*$*$*  RUN  4:  expm  via  eigenvalues  and  eigenvectors  *$*$*$* 

110 

7,  1:  gamma  scaling  I  2:  truncation  I  3:  second  loop  tol  I  4:  unif  ormization 

111 

°l  *****  what  methods  do  we  want  to  run? 

112 

RunData. Method (1) .Use  =  false; 

113 

RunData. Method (2) .Use  =  false; 

114 

RunData. Method (3) .Use  =  true; 

115 

RunData. Method (4) .Use  =  false; 

116 

RunData. Method (5) .Use  =  false; 

117 

7,  *****  what  parameters  needed  for  each  method? 

118 

RunData. Method(l) .Param  =  0; 

119 

RunData. Method (2) .Param  =  0; 

120 

RunData. Method (3) .Param  =  0; 

121 

RunData. Method (4) .Param  =  0; 

122 

RunData. Method (5) .Param  =  0; 

123 

case  5 

124 

7,  *$*$*$*  RUN  4:  Matrix  exponential  via  Taylor  series.  *$*$*$* 

125 

7,  1:  gamma  scaling  I  2:  truncation  I  3:  second  loop  tol  I  4:  unif  ormization 

126 

7,  *****  what  methods  do  we  want  to  run? 

127 

RunData. Method (1) .Use  =  false; 
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128  RunData. Method (2) .Use  =  false; 

129  RunData. Method (3) .Use  =  false; 

130  RunData. Method (4) .Use  =  true; 

131  RunData. Method (5) .Use  =  false; 

132  7,  *****  what  parameters  needed  for  each  method? 

133  RunData. Method (1) .Par am  =  0; 

134  RunData. Method(2) .Param  =  0; 

135  RunData. Method (3) .Param  =  0; 

136  RunData. Method (4) .Param  =  0; 

137  RunData. Method (5) .Param  =  0; 

138  otherwise 

139  7,  *$*$*$*  RUN  6:  expm  via  Pade  approximation  *$*$*$* 

140  7o  1:  gamma  scaling  I  2:  truncation  I  3:  second  loop  tol  I  4:unif  ormization 

141  7o  *****  what  methods  do  we  want  to  run? 

142  RunData. Method (1) .Use  =  false; 

143  RunData. Method (2) .Use  =  false; 

144  RunData. Method (3) .Use  =  false; 

145  RunData. Method (4) .Use  =  false; 

146  RunData. Method (5) .Use  =  true; 

147  7.  *****  what  parameters  needed  for  each  method? 

148  RunData. Method (1) .Param  =  0; 

149  RunData. Method (2) .Param  =  0; 

150  RunData. Method (3) .Param  =  0; 

151  RunData. Method (4) .Param  =  0; 

152  RunData. Method (5) .Param  =  0; 

153  end 

154  disp( [’ Starting  batch  run  for  state  1  int2str (intSt)] ) ; 

155  save  RunData; 

156 

157  7*  call  the  run 

158  tic; 

159  [BestF,BestI ,RimStats ,RunSet]  =  mads_batch; 

160  myoutput  =  toe; 

161 

162  7o  output  results  to  a  text  file 

163  intN  =  [int2str (intSt)  int2str(i)  int2str (rnum)] ; 

164  strFN  =  [ Jmads_output ’  intN]; 

165  7o  output  results  to  a  .mat  file 

166  save (strFN); 

167 

168  myoutput  =  BestF ; 

169 

170  return; 
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Appendix  H.  Post-processing  Code 
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2 
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41 


#/0  Tim  Booher 

7,  2  March  2006  —  This  script  postprocesses  batch  MADS  data  files 

7, 

7,  First,  data  are  conditioned,  then  placed  in  an  excel  file 
clear  all; 

strRpt  =  ’results_f inal . xls ’ ; 

intMethods  =  1:6; 

runCount  =  30 ; 

states  =  [2  5  7  10  20] ; 

stateSize  =  length(states) ; 

workdir  =  pwd; 

lm  =  length (intMethods) ; 

TAU  =  zeros (stateSize, lm) ; 

A  =  zeros (stateSize, lm) ; 

FCALL  =  A; 

C  =  A; 

RTIME  =  A; 

RUNTIMES  =  zeros (1,1m); 

7,  collect  the  data 
for  intCurSt  =  states 

k  =  f ind(states==intCurSt) ; 
f  or  myrun  =  1 : runCount 
mr  =  int2str (myrun) ; 
for  j  =  intMethods 

CurCol  =  f ind(intMethods==j ) ; 
cs  =  int2str (intCurSt) ;  /  current  state 
cm  =  int2str(j);  %  current  method 
mydir  =  [cs  cm  mr] ; 

fname  =  [’ mads _ output’  cs  cm  mr  ’.mat’]; 
disp(fname) ; 
if  (exist (fname) ==2) 
load(fname) ; 

STATE(k) .TAU (myrun, CurCol)  =  BestF.x; 

STATE (k) .A (myrun, CurCol)  =  -BestF.f; 

STATE(k) .C (myrun, CurCol)  =  -BestF.c; 

STATE(k) . FCALL (myrun, CurCol)  =  RunStats .nFunc ; 
STATE(k) . RTIME (myrun, CurCol)  =  RunStats . time ; 

else 

STATE (k) . TAU (myrun , CurCol )  =  NaN ; 

STATE(k) .A (myrun, CurCol)  =  NaN; 
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42 


STATE(k) .C(myrun, CurCol)  =  NaN; 
STATE (k) .FCALL(myrun,CurCol)  =  NaN; 
STATE (k) .RTIME(myrun,CurCol)  =  NaN; 


43 

44 

45  end 

46  end 

47  end 

48  SUMT(k) .TAU.HW  =  1 . 96* (sqrt (var (STATE (k) . TAU) )/sqrt (runCount) ) ; 

49  SUMT(k) .TAU. Mean  =  mean (STATE (k) . TAU) ; 

50  SUMT(k).A.HW  =  1 . 96* (sqrt (var (STATE(k) . A) ) /sqrt (runCount) ) ; 

51  SUMT (k) . A . Mean  =  mean (STATE (k) . A) ; 

52  SUMT(k).C.HW  =  1 . 96* (sqrt (var (STATE(k) . C) ) /sqrt (runCount) ) ; 

53  SUMT (k) . C . Mean  =  mean (STATE (k) . C) ; 

54  SUMT(k) .FCALL  =  mean (STATE (k) . FCALL) ; 

55  SUMT(k) .RTIME.HW  =  1 . 96* (sqrt (var (STATE (k) . RTIME) ) /sqrt (runCount )) ; 

56  SUMT (k) .RTIME. Mean  =  mean (STATE (k) . RTIME) ; 

57  SUMT(k) .TVAR  =  var (STATE (k) .RTIME) ; 

58  RUNTIMES  =  [RUNTIMES;  STATE (k) . RTIME] ; 

59  end 

60 

61  l  output  the  data 

62  AVGM  =  zeros(5,stateSize*lm) ; 

63  RAW  =  zeros (runCount*5, stateSize*lm) ; 

64  EXPM  =  zeros (stateSize*2 , length (intMethods) ) ; 

65  myCol  =  0;  expmrow  =  -1; 

66  for  j  =  states 

67  k  =  f ind(states==j) ; 

68  expmrow  =  expmrow  +  2; 

69  myExpmCol  =  0; 

70  for  i  =  intMethods 

71  curCol  =  f ind(intMethods==i) ; 

72  myCol  =  myCol  +  3; 

73  myExpmCol  =  myExpmCol  +  1 ; 

74  AVGM(l,myCol-2)  =  SUMT(k) .TAU. Mean (curCol) -SUMT (k) .TAU.HW ( curCol ) ; 

75  AVGM (1, myCol- 1)  =  SUMT (k) .TAU. Mean (curCol ) ; 

76  AVGM(1, myCol)  =  SUMT (k) .TAU. Mean (curCol )+SUMT(k) .TAU.HW ( curCol ) ; 

77  AVGM(2,myCol-2)  =  SUMT (k) .A. Mean (curCol) -SUMT (k) .A. HW( curCol) ; 

78  AVGM (2, myCol- 1)  =  SUMT (k) .A. Mean (curCol) ; 

79  AVGM(2, myCol)  =  SUMT (k) .A. Mean (curCol) +SUMT(k) .A. HW( curCol) ; 

80  AVGM(4,myCol-2)  =  SUMT (k) .FCALL (curCol ) ; 

81  AVGM (4, myCol- 1)  =  NaN; 

82  AVGM (4, myCol)  =  NaN; 

83  AVGM(3,myCol-2)  =  SUMT (k) .RTIME. Mean ( curCol ) -SUMT (k) . RTIME.HW (curCol ) ; 

84  AVGM (3, myCol- 1)  =  SUMT (k) . RTIME .Mean (curCol ) ; 
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85  AVGM(3 ,myCol)  =  SUMT(k) .RTIME. Mean ( curCol )+SUMT(k) .RTIME.HW( curCol) 

86  AVGM(5 ,myCol-2)  =  SUMT(k) . TVAR(curCol) ; 

87  AVGM(5,myCol-l)  =  NaN; 

88  AVGM(5 ,myCol)  =  NaN; 

89  EXPM(expmrow,  myExpmCol)  =  SUMT(k) .RTIME.Mean(curCol) ; 

90  EXPM(expmrow+l ,  myExpmCol)  =  sqrt(SUMT(k) .TVAR(curCol) ) ; 

91  stRow  =  1; 

92  for  rc  =  l:runCount 

93  RAW ( stRow, myCol)  =  STATE (k) .TAU(rc,curCol) ; 

94  RAW(stRow+l ,myCol)  =  STATE(k) . A(rc , curCol) ; 

95  RAW(stRow+2 ,myCol)  =  STATE(k) .FCALL(rc, curCol) ; 

96  RAW(stRow+3,myCol)  =  STATE (k) . RTIME(rc , curCol) ; 

97  RAW(stRow+4,myCol)  =  NaN; 

98  stRow  =  stRow  +  5; 

99  end 

100  end 

101  end 

102 

103  xlswrite(strRpt ,  AVGM,  ’summary’,  ’B3’); 

104  xlswrite(strRpt ,  RAW,  ’raw’,  ’B3’); 

105  xlswrite(strRpt ,  EXPM,  ’expm’,  ’B3’); 

106  xlswrite(strRpt ,  RUNTIMES,  ’runtimes’,  ’B3’); 

107 

108  disp( ’done ’ ) ; 
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Appendix  I.  Plot  Generator 

1  function  [A, cr , cd, ci ,ERt ,mym,tau]  =  Controller (intState) 

2  it  ***************  Enumeration  Tau  Plotter  ****************************** 

3  °/0  Tim  Booher  —  timothy.booher@afit.edu 

4  i  051018 

Ej  J  :|c************************************************************************ 

6  °/0intState  =  2;  i  #  of  states  in  random  environment 

7  run_num  =  intState;  i  unique  id  to  save  files  associated  with  this  run 

8  blnUseFlops  =  false; 

9  blnOutputText  =  false;  °/0  create  text  file? 

10  blnProdPlots  =  false;  i  produce  plots 

11  MOC  =2;  i  which  method  are  we  concerned  with?  (in  plots) 

12  inc  =  25;  i  how  many  increments? 

13  mult  =0.3;  °/0  how  far  from  \tau  =  0? 

14  sp  =1;  i  what  range  of  the  overall  \Lambda  do  we  cover? 

15  1  Runlnfo  —  contains  information  on  each  run 

16  it  ********************************************************************* 

17  i  *$*$*$*  RUN  1  *$*$*$* 

18  i  *****  name  for  the  run 

19  Runlnf o(l) .runTitle  =  [num2str (intState)  1  state  with  unif ormization’ ] ; 

20  i  *****  numerical  id  for  the  run 

21  Runlnf o(l) .runID  =  1; 

22  i  1: gamma  scaling  I  2:first  loop  tol  I  3: second  loop  tol  I  4 :unif ormization 

23  i  *****  what  methods  do  we  want  to  run? 


24  Runlnfo(l) .Method(l) .Use  =  true; 

25  Runlnf o(l) .Method(2) .Use  =  false; 

26  Runlnf o(l) .Method(3) .Use  =  false; 

27  Runlnf o(l) .Method(4) .Use  =  false; 

28  i  *****  what  parameters  needed  for  each  method? 

29  Runlnfo(l) .Method(l) .Param  =  le-3; 

30  Runlnf o(l) .Method (2) .Param  =  le-3; 

31  Runlnf o(l) .Method (3) .Param  =  0; 

32  Runlnf o(l) .Method (4) .Param  =  100; 

33 


34  it  ************************************************************************* 

35  intiter  =  0; 

36 

37  if  blnUseFlops 

38  if  (exist (’C:\Program  Files\ICL\WinPAPI\PAPI  MATLAB  Support\f lops . dll 5 ) ~=3) 

39  path( 5C : \Program  Files\ICL\WinPAPI\PAPI  MATLAB  Support\ ’ ,path) ; 

40  disp(’path  added  for  FLOPS’); 

41  else 
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42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 
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63 

64 
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67 
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76 
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84 


disp(’path  addition  not  needed  —  FLOPS  in  use’); 

end 

end 

Jt  3|e3(e3|c3|e3(e3(c3|e3fc3|e3|e3fe3|c3fe3(c3|e3fe3(e3)e3(c3(c3fe3fe3|e3|e3fc3|e3|e3(c3|c3)e3(c9)e3(e3(c3)e3fc3|e3|e3(e3te3(e3|e3|e3(e3(c3|e3fc3(c3ie3f:3|e3(e3|e3|e3|e3(c3)e3(e3(c3fe3f:3(e3|e3(e3|e3(e3(c3|e3fe3(c9|e3(:3|e 

l  create  output  file 
if  (blnOutputText) 

fid  =  f open ([’ output 5  int2str (run_num)  ’ . csv’] , ’wt ’ ) ; 
fid2  =  f open ( [’watch’  int2str (run_num)  ’ . csv’] , ’wt ’ ) ; 
fid3  =  f open( [’p_rpt ’  int2str (run_num)  ’ . csv’] , ’wt ’ ) ; 

end 

3|c3te3(c3|e3fe3(c3(e3fe3|e3|e3fe3|c3fe3(c3|e3(e3(e9|e3)e3(c3|e3(e3|e3|e3fe3|e3fe3(c3|e3(e3(e3|e3(c3(c3|e3fc3|e3fe3fe3|e3fe3|e3|c3(e3(c3(e34e3(c3ie3f:3|e3(e3|e3|e3|e3(c9|e3(e3(c3|e3f:3(c3(e3(e3|e3(e3(c3|e3(e3(c3|e3(:3|e 

7*  initiate  loop  over  all  methods 
for  i  =  1 : length (Runlnfo) 
if  (blnUseFlops) 
f lops(0) ; 

else 

tic; 

end 

^c:f:3|e^c3|c^e^c3|e^e^0|e^c^e3|e^c4;^e^c3|e3|e^c3|c3|e^e3|e^c^c3|e^e^0|e^c3|e^e^c3|e^e^e3|c^c3f;3|e^c^;^e^c3|e^e^e3|c^e3|e3|e^c3f;^e^c3f;^e^c3|c^e3f;3|e^c4:3|e^c^e 

7*  save  Runlnfo  into  .mat  file 
Runlnf o(i) . intState  =  intState; 

RunData  =  Runlnf o(i); 
save  RunData; 

^  ^c3)o|e^c3|c^e^c3|e^e3|o|e^e^e3|e^e^e^e^c3|e3|e^:3|c^e3f;3|e^c^C3|e^c^o|e^c3|e^e^c3|e^e^e3|c^c3f;3|e^c^e^e^c3|e^e^e^c^e3|o|e^c3f;3|e^c3f;^e^c3|c3|e^e3|e^c^e3|e^c^e 

7,  get  inputs 

Param  =  feval( ’availabilityCalc_Param’ ) ; 
setappdata(0, ’PARAM’ , Param) ; 

Lambda  =  sp*Param. lim/min (Param. r) ; 

:j  ^c:f:3|e^c3|c3|e^c3|e^e3|o|e^e^e3|e^c^e^e^c3|c3|e^:3|e^e3f;3|e^c^o|e^c^o|e^c3|e^e^c3|e^e^e3|c^c3f;3|e^c^;^e^c3|e^e^e3|c^e3f;3|e^c3f;3|e^c3f:3|e^c3|c3|e^e3|e^c^e3|e^c^e 

7,  build  tau  vector 

tau  =  linspace (mult*Lambda/inc, Lambda-Lambda/ (inc*mult) ,inc) ; 

7,lt  =  length  (tau)  ;  7«  i  think  this  can  be  ’inc’ 

A  =  []  ;  cr  =  []  ;  cd  =  []  ;  ci  =  []  ;  ERt  =  []  ; 

:j  ^c:f:3|e^c3|c^e^c3|e^e^o|e^e^e3|e^c3f;3|e^c3|e3|e^c3|c3|e3f;3|e^c^o|e^c^o|e^c3|c^e^c3|e^e^e3|c^c3f;3|e^c^;^e^c3|e^e^e^c^e3f;3|e^c3f;^e^c3f;3|e^c3|c3|e^e3|e^c^e3|e^c^e 

7,  Perform  Iterations  over  all  \tau 
mym  =  1 : inc ; 
for  m  =  1 : inc 

if  (Param. lim/ (min (Param. r)*tau(m))  <=  1) 

error(’\tau  is  too  big  relative  to  \Lambda  —  add  increments’); 
break; 

end 

gamma  =  ceil (Param. lim/ (min (Param. r) *tau(m) )) ; 
fprintf(l,’:  7«s  :  ’  ,RimInf o(i)  .runTitle) ; 
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fprintf(l,’  Iteration: \t  °/0g  of  °/0g J  ,m, inc)  ; 

fprintf (1,’  \t  tau:  7«g  \t  gamma:  70g\n’  ,tau(m) ,  gamma); 

[A(m),  cr(m),  cd(m) ,  ci(m),  ERt(m)]  =  availabilityCalc (tau(m) ) ; 

end 

^  ^c:f:3|e^c3|c^e^c3|e^e^o|e^e^e3|e^c^e^e^c3|o|e^c3|e3|e3f;3|e^c^e3|e^e3f;3|e^c3|c^e^c3|e3|c^e3|c^c^e^e^c^;3|e^c3|e^e^e^c^e3f;3|e^c^e^e^c3f;^e^c3|e^e^e3|e^e^e3|e^c^e 

7*  Output  Performance  Metrics  to  Systerm 
if  (~blnUseFlops) 
ops(i)  =  toe; 

else 

[ops(i),  mflops(i)]  =  flops; 

end 

^  ^c:f;3|e^c3|c^e^c3|e^e^o|e^c^e3|e^c^;^e^c3|o|e^c3|e^e^e3|e^c^e3|c^c3f;3|c^c3|c3|e^c3|e3|c^e^c^c3|o|e^c^e^e^c3|e3|e^e3|c^e3f;3|e^c^e^e^c3f;3|e^c3|e^e^e3|e3|e^e3|e^c^e 

7*  calculate  performance  on  any  iteration  other  than  1  (baseline) 
for  j  =  2:i 

pim  =  100*((ops(l)-ops(j))/ops(l)) ; 
mse  =  mean( (AC{j}-AC{l}) . ~2) ; 
if  (blnOutputText) 

fprintf  (fid3, ’7«g,  °/0g,  ’  ,  run_num,  pim); 
fprintf  (fid3, ’7«g,  7og’,  mse,  intiter); 

end 

end 

^  ^c:f;3|e^c3|c^e^c3|e^e^o|e^c^e3|e^c^;^e^c3|o|e^:3|e^e3f;3|e^c^e3|e^e^o|c^c3|c^e^c3|e3|e^e^c^e^o|e^c^e^e^c3|e3|e^e3|c^e3f;3|e^c^e^e^c3f;3|e^c3|e^e^e3|e3|e^e3|c^c^c 

7,  Output  Data  to  text  file  for  post-processing  —  used  on  remote  unix 
7,  system  (HPC) 
if  (blnOutputText) 
for  b  =  l:inc 

fprintf  (fid,  J7og}  ,tau(b)  )  ; 
fprintf  (fid,  ’ ,  7og’,A(b)); 
fprintf  (fid,  ’ ,  7»g’,C(b)); 

intM  =  [Runlnf o(i) .Method(l) .Use  Runlnfo(i) .Method(2) .Use] ; 
intM  =  [intM  Runlnf o(i) .Method(3) .Use  Runlnf o(i) .Method (4) .Use] 
fprintf  (f  id,  ’ ,  7»s  ’  ,mat2str  (intM)  )  ; 
fprintf  (f  id, } ,  7og’  ,run_num)  ; 
fprintf (f id , J \n ’ ) ; 

end 

save( [’ output int2str(run_num)] ,  5 A’ ,  ’  C’,  ’tau’); 

end 

end 

3|e^c^e3|e^c3f:3|e^c3|c3|e^c3|e^e3f;3|e^c3|e3|e^c^e^e^c3|e3|e^:3|c^e3f;3|e^c^o|e^c^o|e^c3|e^e^c3|e^e^e3|c^c3f;3|e^c^;^e^c3|e^e^e3|c^e3f;3|e^c3f;^e^c3f;^e^c3|c3|e^e3|e^e^e3|e^c^e 

7,  Display  final  results 

3|e3)e3|e3|e3fe3(e3|e3fe3|c3fe3fe3|e3(e3|c3|e3fe3(e3|e3)e3(c3fe3fe3|e3|e3fc3|c3fe3(c3|e3(e3|c3|e3fe3(c3|e3fc3|c3fe3fe3|c3(e3|e3|e3fe3(c3|e3fe3(c3fe3fe3)e3(e3|e3(e3fe3(c9|e3(:3(c3fe3fe3(c3fe3fe3|e3fe3(e9|e3(e3|e3|e3)e3|e 

7,  Produce  Plots 
if  (blnProdPlots) 
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figure;  hold  on; 
myLineSpec{l}  =  ’-k’; 
myLineSpec{2}  =  ’-or’; 
myLineSpec{3}  =  ’-+g’; 
myLineSpec{4}  =  * — xc’ ; 
myLineSpec{5}  =  ’ — sm’ ; 
myLineSpec{6}  =  ’ — db’ ; 
la  =  length (AC) ; 
if  (la  >  6) 

warning ( ’More  methods  used  than  linespecs  available.’); 

end 

for  j  =  1:1a 

plotyy(tau,-AC{j},tau,CC{j}) ; 

end 

ylabel( ’Availability (\tau) ’ ) ; 
xlabel(’\tau’) ; 

*/,  legend(plts , legStr) ; 

title( [int2str (intState)  ’-state  case  with  ’... 

int2str(inc)  ’  increments’]); 
grid  on; 

myx  =  get (gca, ’xlim’ ) ;  myy  =  get (gca, ’ylim’ ) ; 
myx  =  myx(l)+(myx(2)-myx(l))*(l/2) ; 
myy  =  myy ( 1) + (myy (2) -myy (1))* (1/2) ; 

:j  ^c:f:3|e^c3|c^e^c3|e^e3|o|e^c^e3|e^c3f;^e^c3|e3|e^c3|c3|e3f;3|e^c^o|e^c^o|e^c3|e^e^c3|e^e^e3|c^c3f;3|e^c^e^e^c3|e^e^c^c3|e^e3|e^c3f;^e^c3f;^e^c3|c3|e^e3|e^c^o|e^c 

°/*  calculate  display  improvement 

^:3f;3|e^c3|e3|e^c3|c^e3f;3|e^e^e3|c^c^o|e^c3|e3|e^:3|e3|e^e3|e^c^e3|c^c3|o|e^:3|ofe:4eHe:ic:le3le:4C3ii:3le:4c:(C3ic:4e3le,(e:le3le9fe3f:He3lC3(C9fe:4e3f;3le3|eHo|eHi:3le:le3ii:3k:4c 

for  j  =  2:i 

text (myx, myy , [’MSE  Method  ’  num2str(j)  ’  =  ’  num2str (mse)] , . . . 
’BackgroundColor ’ , ’white’) ; 

text (myx, myy, [’Improvement  Method  ’  num2str(j)  ’  =  ’  num2str (pim)] , . . . 

’BackgroundColor’ , ’white’) ; 
text (myx, myy, [’Method  ’  num2str(j)  ’  Prm  =  ’  ... 

num2str (Runlnf o(j) . Method (MOC) .Param)] , ’BackgroundColor’ , ’white’) ; 

end 

end 

if  (blnOutputText) 

fclose(fid);  f close(f id2) ;  f close(f id3) ; 

end 

return; 
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Appendix  J.  Perl  Batch  Controller 


1 

#! /bin/perl  -w 

2 

use  Cwd;  use  File:: Copy; 

3 

my  $stdir  =  get cwd () ;  #  get  current  directory 

4 

print  "starting  dir  is  $stdir\n" ; 

5 

my  $nr  =  1;  my  $i  =  0;  #  initialize  index 

6 

foreach  $st  (qw/2  5  7  10  20/)  {  #  iterate  over  all  states 

7 

@j  =  qw/5  10/;  #  array  of  times  for  each  state 

8 

for  ($k=l;  $k  <=  6;  $k++)  {  #6  different  methods 

9 

for  ($r=l;  $r  <=  $nr;  $r++)  {  #  iterate  over  each  run 

10 

$r  =  sprintf  ("°/002d"  ,  $r)  ;  #  format  the  last  two  digits 

11 

$rnum=$st . $k. $r ;  #  unique  id  for  each  run 

12 

print  "submitting  $i:$st  \t  $rnum  \t  k:  ($k)  \t  $j[$i]\n"; 

13 

$myscript  =  <<"ENDSCRIPT" ; 

14 

# ! /bin/sh 

15 

#BSUB  -q  regular 

16 

#BSUB  -n  1 

17 

#BSUB  -W  $j  [$i] : 00 

18 

#BSUB  -a  SMP 

19 

#BSUB  -J  b_$rnum 

20 

#BSUB  -o  matlab_run_$rnum. out 

21 

#BSUB  -e  matlab_run_$rnum. out 

22 

#BSUB  -P  WPBAFIT025047ACM 

23 

mat lab  -nojvm  -nodesktop  -nosplash  >  matlab_run_$rnum . out  <<  EOL 

24 

CallMadsBatch($st , $k, $r , $stdir) ; 

25 

EOL 

26 

ENDSCRIPT 

27 

mkdir ($rnum,0777) ;  #  create  a  directory  with  full  permissions 

28 

system  "cp  *.m  $rnum/"; 

29 

chdir  $rnum;  #  place  all  files  in  that  directory 

30 

$myfilename  =  "booher_job_$rnum. sh" ; 

31 

open  DUMMYFILE,  ">$myf ilename" ;  print  DUMMYFILE  $myscript; 

32 

close  DUMMYFILE; 

33 

system  "bsub  <  $myf ilename" ;  #  submit  the  script  to  the  lsf  engine 

34 

chdir  $stdir; 

35 

#print  $rnum. "\n" ; 

36 

> 

37 

> 

38 

$i+=l; 

39 

> 

j-i 
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